text:
# An interest rates cluster analysis

## 1 Introduction

Since long time financial data have been widely studied by economists, mathematicians and, more recently, by physicists . The variations of these financial time series can be seen as stochastic processes where a set of financial quantities is varying in time as consequence of underlying economic changes. The present availability of enormous sets of financial data allows to get increasingly important insights on the complex behavior of these systems starting from empirical studies. These investigations are leading to more and more accurate results on risk assessment and search for market imperfections. One of the important points in these analyses is to individuate similarities and specificities among the analyzed financial time series. This search has been widely exploited for stocks price changes whereas interest rates have been less investigated . For several economic reasons, interest rates and bonds have very similar statistical behavior in time or, in other words, they are all highly correlated . Their multivariate dynamics has been studied with a correlation-based clustering procedure in a set of US treasury securities where an underlying hierarchical structure has been detected . Here we investigated a partially different and less homogeneous set to evaluate the degree of hierarchal organization among different time series in a diversified group of bonds.

## 2 An empirical analysis on interest rates

We investigate weekly data for 34 selected interest rate time series recorded in the Federal Reserve (FR) Statistical Release database . In the following we will indicate these time series with the symbol $`f_i(t)`$, where $`t`$ is the current date and $`i`$ is a number which labels the different time series (see Tab.1).

The different interest rate time series analyzed are: The Federal funds rate (FED); State $`\&`$ local bonds (SLB); Commercial Paper (CP); Finance Paper placed directly (FP); Bankers acceptances (BA); The rate on certificates of deposit (CD); (Note that in these cases the numbers $`1`$, $`3`$ and $`6`$ stand for maturity dates of $`1`$, $`3`$ and $`6`$ months.); The yields on Treasury securities at ‘constant maturity’ (TC) (in particular the TC at 3 and 6 months (TC3M, TC6M) and 1, 2, 3, 5, 7, 10, and 30 years (TC1Y-TC30Y) maturities; The Treasury bill rates (TBA) with maturities of 3 and 6 months (TBA3M, TBA6M, TBS3M, TBS6M) and 1 year (TBS1Y); The Treasury long-term bond yield (TC10P); The Eurodollar interbank interest rates (ED) with maturity dates 1 month, 3 months and 6 months (ED1M, ED3M, ED6M), respectively; The Corporate bonds Moody’s seasoned rates (AAA, BAA) and The Conventional mortgages rates (CM). Their characteristics can be found in . Unless differently stated, we report weekly data obtained from unweighted averages of daily data ending on Friday.

## 3 Fluctuations

The interest rate time series, $`f_i(t)`$ v.s. $`t`$ are shown in Fig.1, where their average $`\overline{f}(t)=_if_i(t)/34`$ is also shown. It is evident from Fig.1 that all these data follow very similar trends in time and they lay in a narrow band around $`\overline{f}(t)`$.

The interest rate fluctuations are analyzed by studying the changes in their values from one week to the following week: $`\mathrm{\Delta }f_i(t)=f_i(t+\mathrm{\Delta }t)f_i(t)`$ where $`\mathrm{\Delta }t=1`$ week. The quantities $`\mathrm{\Delta }f_i(t)`$ show stochastic fluctuations around the zero with similar behaviors for all the interest rates. These fluctuations are analyzed in the ‘tails’ region by computing the cumulative distributions $`\mathrm{\Psi }_\pm (\pm x)`$ ($`x>0`$), a quantity which tells us the probability to find a weekly change which is *larger than* $`x`$ ($`+`$), or *smaller than* $`x`$ ($``$). It is defined as: $`\mathrm{\Psi }_+(x)=1_{\mathrm{}}^xp(\xi )𝑑\xi `$ and $`\mathrm{\Psi }_{}(x)=_{\mathrm{}}^xp(\xi )𝑑\xi `$ with $`p(\xi )`$ being the probability density distributions of $`\mathrm{\Delta }f_i(t)`$ (see Fig.1 (Right)). These distributions are highly leptokurtic and are characterized by non-Gaussian profiles. The standard deviation of $`\mathrm{\Delta }f_i(t)`$ is defined as: $`\sigma _i=\sqrt{\frac{1}{T_2T_1}_{t=T_1}^{T_2}(\mathrm{\Delta }f_i(t)\mathrm{\Delta }f)^2}`$ where $`T_1`$ and $`T_2`$ delimit the range of $`t`$, and $`\mathrm{\Delta }f`$ is the average over time of $`\mathrm{\Delta }f_i(t)`$ (which tends to zero for $`T_2T_1\mathrm{}`$). We compute the standard deviations of $`\mathrm{\Delta }f_i`$ for each interest rate series for the whole period 1982-1997 (see Tab.1) and for each year (see Fig.2 (Left)). We can observe an overall decreasing trend of $`\sigma _i`$ in the time period 1982-1997 with similar fluctuations for all the interest rates series except for FED.

## 4 Cluster analysis and discussion

To understand the geometrical and topological structure of the correlation coefficients, we use the metric distance $`d_{i,j}`$ between the series $`\mathrm{\Delta }f_i`$ and $`\mathrm{\Delta }f_j`$ which is defined in  and used for financial time series in : $`d_{i,j}=\sqrt{2(1c_{i,j})}`$ with $`c_{i,j}`$ the correlations among the $`i,j`$ interest rates weekly changes:

$$c_{i,j}=\frac{\mathrm{\Delta }f_i\mathrm{\Delta }f_j\mathrm{\Delta }f_i\mathrm{\Delta }f_j}{\sigma _i\sigma _j},$$
(1)
where the symbol $`\mathrm{}`$ denotes a time average performed over the investigated time period. The correlation coefficients are computed between all the pairs of indices labeling our interest series. Therefore we have a $`34\times 34`$ symmetric matrix with $`c_{i,i}=1`$ on the diagonal. By definition, $`c_{i,j}`$ is equal to zero if the interest rates series $`i`$ and $`j`$ are totally uncorrelated, whereas $`c_{i,j}=\pm 1`$ in the case of perfect correlation/anti-correlation. Therefore $`d_{i,j}`$ can vary between 0 to 2. We determine an ultra-metric distance $`\widehat{d}_{i,j}`$ which satisfies the first two properties of the metric distance and replaces the triangular inequality with the stronger condition: $`\widehat{d}_{i,j}max[\widehat{d}_{i,k},\widehat{d}_{k,j}]`$, called ‘ultra-metric inequality’. Once the metric distance $`d_{i,j}`$ is used, one can introduce several ultra-metric distances. Mantegna et. al have used the ‘subdominant ultra-metric’, obtained by calculating the minimum spanning tree connecting several financial time series . Here we consider a different ultra-metric space that emphasizes the cluster-structure of the data. In our case, a “cluster” is a set of elements with relative distances $`d_{i,j}`$ which are smaller than a given threshold distance $`\overline{\delta }`$, whereas disjoined clusters have some elements which are at distances larger than $`\overline{\delta }`$. We define the *ultra-metric distance* $`\widehat{d}_{i,j}`$ between two distinct elements $`i,j`$ belonging to two different clusters as the maximum metric distance between all the couples of elements in the two clusters . The linkage procedure yields to a hierarchical graph as shown in Fig.2 (Right), which refers to the 34 time series in the whole period 1982-1997. The clustering process starts with a nucleation between TC1Y and TBS1Y at the ultra-metric distance $`\delta _n=0.087`$. The clustering ends when all the series merge in a unique large cluster at the ultra-metric distance $`\delta _p=1.3083`$. At this distance all the interest rates with maturity dates smaller or equal than 6 months (already merged with the FED at $`\widehat{d}=1.18`$) make a single cluster with another cluster composed of interest rates with maturity dates larger or equal than 1 year. The same analysis performed on each year, gives comparable $`\delta _n`$ and $`\delta _p`$ values which are plotted in Fig.3 (Left). Let now consider the intermediate region by analyzing the cluster evolution at the threshold distance $`\overline{\delta }=1/\sqrt{2}=0.707\mathrm{}`$, which is half the way between completely uncorrelated series ($`c_{i,j}=0`$ and $`d_{i,j}=\sqrt{2}`$) and completely correlated ones ($`c_{i,j}=1`$ and $`d_{i,j}=0`$). At this threshold distance, the cluster analysis on the whole data set (1982-1997) leads to 6 clusters and 3 isolated elements, as one can see from Fig.2 (Right).

The corresponding interest rates associated with these clusters are summarized in first column (1982-1997) of Fig.3 (Right) where the ultra-metric distances at which the nucleation process starts and ends for each cluster are also indicated. As one can see, the empirical analysis allow us to distinguish several different clusters that gather together meaningful quantities:

* all the interest rates with maturities equal to 1 months;
* all the interest rates with maturities 3 and 6 months with distinctions for the Treasury securities at ‘constant maturity’ (TC), Treasury bill rates (TBA) and Treasury bill secondary market rates (TBS);
* all the interest rates with maturities between 1 and 3 years;
* all the interest rates with maturities larger than 3 years.

Fig.4 reports the plot for the interest rates time series $`f_i(t)`$ grouped into the different sets retrieved from the cluster analysis described above. For some of the series the data-collapse is impressive, indicating that the correlations inside the clusters are strong in any part of the analyzed period. It is therefore interesting to investigate whether a cluster structure, similar to the one obtained for the period 1982-1997, could be retrieved from an analysis on shorter periods. This is of course a delicate point since the fragmentation of the data sets will increase the fluctuations due to the noise. We choose to divide the whole period 1982-1997 in four smaller periods of 4 years. The results are reported in Fig.3 (Right) where it is evident how the cluster structure is mostly conserved (and partially modified) in this 4-years period analysis. In conclusion, from the analysis of different kinds of interest rates in money and capital markets, referring to government, private, industries securities and commitments, we have shown how the used clustering linkage procedure is useful to detect differences and analogies among these tangled correlated data.


format: text

text:
# The influence of complex material coverings on the bandwidth of antennas

## 1 Introduction

Properties of antennas depend on materials that may cover metal radiating parts or fill the antenna volume. For instance, it is known that the bandwidth of microstrip antennas can be improved by using magnetic substrates . In the recent literature there is renewed interest to the question of improving antenna performance by loading the antenna with complex materials. In particular, it is of interest if new artificial materials with negative material parameters  can be useful in the antenna design. This problem has been analyzed in a recent conference presentation  for a point dipole antenna surrounded by a shell of a material with negative parameters, and some very promising results reported. However, the key issue of the reactive energy storage and the antenna quality factor calls for a more careful examination, because of inevitable frequency dispersion of passive artificial materials.

In this paper we analyze in detail a very simple radiating system: an infinite current line in a material cylindrical shell. Assuming the diameter to be small compared with the wavelength in the material, the field equations can be solved analytically leading to simple and physically clear formulas. The quality factor can be then defined and calculated for arbitrary dispersive and also for specific lossy coverings. The results of the analysis lead to some general conclusions regarding the effect of thin material coverings on the antenna performance.

## 2 Fields of a thin radiating material-covered cylinder

We consider the radiation from a line electric or magnetic current embedded into a material cylinder (Figure 1). The problem originates from a simplified model of a thin-wire antenna covered by a thin layer of some material. The antenna bandwidth is investigated by means of the effective impedance per unit length of the antenna wire. The real part of the introduced effective impedance corresponds to the radiated power (radiation resistance), and the imaginary part gives information about the energy stored in the antenna near field. A thin material layer covering the radiating current affects the stored energy and, if the material is lossy, the radiated power. Our aim is to find how this changes the antenna bandwidth.

Let us consider an infinitely thin line of time-harmonic electric current $`I`$. We assume the current to be uniform in its phase and its amplitude along the line. We model a thin radiating wire of radius $`r_0`$ as such current line sitting at the wire axis. This method is similar to taking into account only the local (singular) part of the thin-wire kernel (e.g., ). We denote the covering cylinder outer radius as $`a`$, the frequency of operation as $`\omega `$, and the material parameters of the outer space and the covering as $`ϵ_0`$, $`\mu _0`$, and $`ϵ`$, $`\mu `$, respectively. Hence, for the wave number in the outer space we have $`k_0=\omega \sqrt{ϵ_0\mu _0}`$, and for that in the material we have $`k=\omega \sqrt{ϵ\mu }`$. The branch of the square root in the last relation is chosen so that $`\mathrm{Im}\{k\}0`$ (the time dependence is of form $`e^{+j\omega t}`$).

Expressing the line source fields inside and outside the cover and using the assumed cylindrical symmetry we write for the electric field:

$$\{\begin{array}{cc}E_z=AH_0^{(2)}(kr)+BH_0^{(1)}(kr),r_0ra,\hfill & \\ E_z=CH_0^{(2)}(k_0r),ra,\hfill & \end{array}$$
(1)
and for the magnetic field:

$$H_\phi =\frac{j}{\omega \mu }\frac{E_z}{r}=\{\begin{array}{cc}\frac{jk}{\omega \mu }\left(AH_1^{(2)}(kr)+BH_1^{(1)}(kr)\right),r_0ra\hfill & \\ \frac{jk_0}{\omega \mu _0}CH_1^{(2)}(k_0r),ra\hfill & \end{array}.$$
(2)
Here the $`z`$-axis is directed along the wire and $`r`$ is the radial distance counted from the axis of the wire \[$`(r,\phi ,z)`$ form a usual set of cylindrical coordinates\].

On the surface $`r=a`$ we have the continuity boundary conditions for $`E_z`$ and $`H_\phi `$ field components:

$$AH_0^{(2)}(ka)+BH_0^{(1)}(ka)=CH_0^{(2)}(k_0a),$$
(3)
$$\frac{\mu _0k}{\mu k_0}\left(AH_1^{(2)}(ka)+BH_1^{(1)}(ka)\right)=CH_1^{(2)}(k_0a).$$
(4)
In addition, we use a relation between the wire current $`I`$ and the magnetic field $`H_\phi `$ at the wire surface $`r=r_0`$:

$$2\pi r_0H_\phi =I,$$
(5)
obtaining

$$\frac{2\pi jkr_0}{\omega \mu }\left(AH_1^{(2)}(kr_0)+BH_1^{(1)}(kr_0)\right)=I.$$
(6)
Eqs. (3), (4), and (6) form a complete system of equations. Instead of solving this system directly in terms of Hankel functions we make a simplification assuming that the following thin-wire condition is satisfied: $`kr_01`$. Under this assumption

$$H_1^{(1,2)}(z)\frac{2j}{\pi z},$$
(7)
and Eq. (6) becomes

$$B=A+\frac{\omega \mu }{4}I.$$
(8)
We substitute the last relation into Eqs. (3) and (4) to obtain

$$2J_0(ka)AH_0^{(2)}(k_0a)C=\frac{\omega \mu }{4}IH_0^{(1)}(ka),$$
(9)
$$2J_1(ka)A\frac{\mu k_0}{\mu _0k}H_1^{(2)}(k_0a)C=\frac{\omega \mu }{4}IH_1^{(1)}(ka).$$
(10)
Finally, we get the solution in the following form:

$$A=\frac{\omega \mu I}{8}\frac{H_1^{(1)}(ka)H_0^{(2)}(k_0a){\displaystyle \frac{\mu k_0}{\mu _0k}}H_0^{(1)}(ka)H_1^{(2)}(k_0a)}{J_1(ka)H_0^{(2)}(k_0a){\displaystyle \frac{\mu k_0}{\mu _0k}}J_0(ka)H_1^{(2)}(k_0a)},$$
(11)
$$C=\frac{j\omega \mu I}{4}\frac{J_0(ka)Y_1(ka)J_1(ka)Y_0(ka)}{J_1(ka)H_0^{(2)}(k_0a){\displaystyle \frac{\mu k_0}{\mu _0k}}J_0(ka)H_1^{(2)}(k_0a)}.$$
(12)
The last unknown $`B`$ can be found from Eq. (8).

The obtained solution is quite difficult to analyze in the general case. However, since we are interested in the material effects in the case of thin coverings, when the antenna cross section remains small in terms of the wavelength, we can assume that $`k_0a1`$, $`|k|a1`$. Then, using the known Bessel function asymptotics for small arguments

$$J_0(z)1,J_1(z)\frac{z}{2},Y_0(z)\frac{2}{\pi }\mathrm{log}\frac{\gamma z}{2},Y_1(z)\frac{2}{\pi z},$$
(13)
(here $`\gamma 1.781`$ is the Euler constant), the relations for amplitudes $`A`$, $`B`$ and $`C`$ reduce to

$$A=\frac{\omega I}{8}\left\{\mu _0+\mu +\frac{2j}{\pi }\left[\mu \mathrm{log}\frac{\gamma ka}{2}\mu _0\mathrm{log}\frac{\gamma k_0a}{2}\right]\right\},$$
(14)
$$B=\frac{\omega I}{8}\left\{\mu _0\mu +\frac{2j}{\pi }\left[\mu \mathrm{log}\frac{\gamma ka}{2}\mu _0\mathrm{log}\frac{\gamma k_0a}{2}\right]\right\},$$
(15)
$$C=\frac{\omega \mu _0I}{4}.$$
(16)
In the next section we will use Eqs. (14)–(16) to find the energy stored in the source reactive field.

## 3 Antenna quality factor

The antenna bandwidth is related with the antenna quality factor. The well-know definition for the quality factor of a lossy (in the antenna case losses are mostly due to radiation) resonator is

$$Q=\frac{\omega W}{P},$$
(17)
where $`W`$ is the average reactive energy stored in the resonator, and $`P`$ is the total dissipated power (for an antenna this is a sum of the radiated power and the thermal loss power). Although $`Q`$ defined by this formula is mostly used for resonant systems, it can be also used as a relative measure of ability of an aperiodic system to store reactive energy.

Usually, the quality factor of antennas is evaluated integrating the reactive field energy density over the whole space and integrating the Poynting vector in the far zone. In our present case, where both the stored energy and the radiated energy are due to the current in the antenna wire, we can use a simpler method introducing the effective impedance per unit length of the wire as follows:

$$Z=\frac{E_z(r_0)}{I},$$
(18)
where $`E_z(r_0)`$ means the total longitudinal electric field component created by the considered line current at the surface $`r=r_0`$. The minus sign is because we suppose $`I`$ to be given (external) current. This approach is sometimes called the induced electromotive force method. The field component of interest is given by

$$E_z(r_0)=AH_0^{(2)}(kr_0)+BH_0^{(1)}(kr_0).$$
(19)
Using the thin-wire approximation we express this as

$$E_z(r_0)=A\left(1\frac{2j}{\pi }\mathrm{log}\frac{\gamma kr_0}{2}\right)+B\left(1+\frac{2j}{\pi }\mathrm{log}\frac{\gamma kr_0}{2}\right).$$
(20)
Substituting $`A`$ and $`B`$ from Eqs. (14), (15), and using Eq. (18) we get

$$Z=\frac{\omega \mu _0}{4}+j\frac{\omega }{2\pi }\left[\mu \mathrm{log}\frac{a}{r_0}+\mu _0\mathrm{log}\frac{4}{\gamma k_0a}\right]=\frac{\omega \mu _0}{4}+j\frac{\omega }{2\pi }\left[\mu _0\mathrm{log}\frac{4}{\gamma k_0r_0}+(\mu \mu _0)\mathrm{log}\frac{a}{r_0}\right].$$
(21)
This formula has a clear physical meaning. $`\frac{\omega \mu _0}{4}`$ is the radiation resistance per unit length. The first member in the brackets is the reactance per unit length in the absence of covering. The last member, proportional to $`\mu \mu _0`$, measures the influence of the covering on the reactive input impedance. Note that under our assumptions the radiation resistance does not depend on the parameters of the thin covering cylinder.

The obtained impedance is a complex number: $`Z(\omega )=R(\omega )+jX(\omega )`$, where

$$R(\omega )=\frac{\omega \mu _0}{4}+\frac{\omega \mu ^{\prime \prime }(\omega )}{2\pi }\mathrm{log}\frac{a}{r_0},$$
(22)
$$X(\omega )=\frac{\omega }{2\pi }\left[\mu _0\mathrm{log}\frac{4}{\gamma k_0r_0}+(\mu ^{}(\omega )\mu _0)\mathrm{log}\frac{a}{r_0}\right].$$
(23)
Here we have assumed loss only in the covering material: $`\mu =\mu ^{}j\mu ^{\prime \prime }`$.

### 3.1 Negligible losses in the covering material

Let us first assume that the losses in the material cover can be neglected. In this case we can find the quality factor in terms of the derivative of $`X(\omega )`$ with respect to the frequency. Here, we follow the approach used in . Instead of considering steady-state harmonic oscillations, we study a transient regime characterized by an exponentially growing (or decaying) amplitude of harmonic oscillations. For this regime we introduce a complex frequency:

$$\mathrm{\Omega }=\omega j\alpha .$$
(24)
Positive values of $`\alpha `$ result in growing oscillations, negative values correspond to decaying ones. We assume that $`|\alpha |\omega `$. This implies a very slow rate of the amplitude change, hence, for the reactive energy $`W`$ (averaged over a period of oscillations) stored in the circuit we can write

$$\frac{W}{t}2\alpha W.$$
(25)
The coefficient 2 is there because the energy is proportional to the square of the oscillation amplitude. The circuit impedance as a function of the introduced complex frequency becomes

$$Z(\mathrm{\Omega })=R(\mathrm{\Omega })+jX(\mathrm{\Omega })R(\omega )+\alpha \frac{X(\omega )}{\omega }+j\left[X(\omega )\alpha \frac{R(\omega )}{\omega }\right].$$
(26)
We see that exponential growth of oscillation amplitude results in an additional resistance as well as an additional reactance.

In the case of negligible dissipation loss which we consider at this stage it is possible to split the total active power (proportional to the real part of the effective impedance) into parts corresponding to radiation and increase of the stored reactive energy. Indeed, the radiated power can be calculated knowing the far-zone fields. If there is no dissipation loss in the antenna, then the way how to find the energy storage power is straightforward. A complexity comes when the loss and storage processes are combined (see the next section).

Comparing Eqs. (16) and (22) we see that for our case the radiation resistance (per unit length) is

$$R_{\mathrm{rad}}=\frac{\omega \mu _0}{4}.$$
(27)
Then the energy storage power is given by (we use the effective values of the fields and the current)

$$\frac{W}{t}=\left(\mathrm{Re}\{Z(\mathrm{\Omega })\}R_{\mathrm{rad}}\right)|I|^2=\alpha \frac{X(\omega )}{\omega }|I|^2.$$
(28)
Comparing this with Eq. (25) we find the averaged stored reactive energy (per unit length of the wire)

$$W=\frac{1}{2}\frac{X(\omega )}{\omega }|I|^2,$$
(29)
and obtain for the quality factor:

$$Q=\frac{\omega W}{P}=\frac{\omega }{2R_{\mathrm{rad}}}\frac{X(\omega )}{\omega },$$
(30)
because the total loss power is the radiated power $`P=R_{\mathrm{rad}}|I|^2`$. Differentiating $`X(\omega )`$ given by (23), we get

$$\frac{X(\omega )}{\omega }=\frac{1}{2\pi }\left\{\frac{(\omega \mu )}{\omega }\mathrm{log}\frac{a}{r_0}+\mu _0\left[\mathrm{log}\frac{4}{\gamma k_0a}1\right]\right\}.$$
(31)
We substitute the last result into Eq. (30) and do simple algebra to get finally

$$Q=\frac{1}{\pi }\left\{\frac{(\omega \mu _r)}{\omega }\mathrm{log}\frac{a}{r_0}+\mathrm{log}\frac{4}{\gamma k_0a}1\right\},$$
(32)
where we have introduced $`\mu _r=\mu /\mu _0`$ as the relative permeability of the covering material.

### 3.2 Lossy and dispersive covering cylinder

If the covering material is dispersive and the dissipation losses cannot be neglected, the quality factor cannot be determined from the analysis of the antenna impedance only. The reason is that the active part of the impedance determines the total active power delivered to the antenna, but we need to know what part of this power is used to increase the stored energy and what part is dissipated into heat. The quality factor can be, however, found, if we actually know the internal structure of the analyzed system. In our particular case, the result will depend on the structure and design of the loading material.

Let us assume that the magnetic material that fills the covering cylinder can be modeled by a Lorentzian permeability:

$$\mu =\mu _0\left(1+\frac{A\omega ^2}{\omega _0^2\omega ^2+j\omega \mathrm{\Gamma }}\right),$$
(33)
where the magnitude factor $`A`$ and the loss factor $`\mathrm{\Gamma }`$ do not depend on the frequency. This material can be realized, for example, as a dense array of small split-ring resonators made of metal wires or strips. The total impedance of the covered antenna (21) we write as a sum of the antenna impedance without material cover and an additional impedance due to the material cylinder:

$$Z=Z_{\mathrm{wire}}+Z_{\mathrm{medium}}$$
(34)
with

$$Z_{\mathrm{wire}}=\frac{\omega \mu _0}{4}+j\frac{\omega \mu _0}{2\pi }\mathrm{log}\frac{4}{\gamma k_0r_0},Z_{\mathrm{medium}}=j\frac{\omega }{2\pi }(\mu \mu _0)\mathrm{log}\frac{a}{r_0}.$$
(35)
If the permeability obeys (33), to model $`Z_{\mathrm{medium}}`$ we can introduce a magnetic-coupled circuit shown in Figure 2. Indeed, in terms of the introduced effective parameters, the additional impedance associated with the material can be written as

$$Z_{\mathrm{material}}=\frac{j\omega ^3(M^2/L)}{\omega _0^2\omega ^2+j\omega (R/L)}$$
(36)
(here $`\omega _0=1/\sqrt{LC}`$) which has the same form as the impedance following from (35) and (33). Given the total current through the antenna wire $`I`$, we can now find the currents and voltages on the reactive elements of the equivalent circuit and calculate the averaged stored energy as

$$W_{\mathrm{medium}}=L\frac{|I_L|^2}{2}+C\frac{|U_C|^2}{2}=\frac{1}{2}\left(L+\frac{1}{\omega ^2C}\right)|I_L|^2.$$
(37)
The result is

$$W_{\mathrm{medium}}=\frac{\omega ^2M^2C(1+\omega ^2LC)}{(1\omega ^2LC)^2+\omega ^2R^2C^2}\frac{|I|^2}{2}.$$
(38)
At the resonant frequency the last expression becomes

$$W_{\mathrm{medium}}(\omega _0)=\frac{M^2|I|^2}{R^2C}.$$
(39)
We observe that if there is no coupling, the additional energy stored in the medium is zero, as it should be. The dissipated power is

$$P_{\mathrm{loss}}=R|I_L|^2=\frac{\omega ^4M^2C^2R|I|^2}{(1\omega ^2LC)^2+\omega ^2R^2C^2}$$
(40)
At the resonant frequency the dissipation is given by

$$P_{\mathrm{loss}}(\omega _0)=\frac{\omega _0^2M^2|I|^2}{R}.$$
(41)
The radiation efficiency of the considered system at the resonant frequency can be found as

$$\eta =\frac{P_{\mathrm{rad}}}{P_{\mathrm{rad}}+P_{\mathrm{loss}}}=\frac{1}{1+{\displaystyle \frac{4\omega _0M^2}{\mu _0R}}},$$
(42)
and the total antenna quality factor is

$$Q_{\mathrm{total}}=\frac{\omega (W_{\mathrm{wire}}+W_{\mathrm{medium}})}{P_{\mathrm{rad}}+P_{\mathrm{loss}}}=\eta \frac{\omega (W_{\mathrm{wire}}+W_{\mathrm{medium}})}{P_{\mathrm{rad}}}.$$
(43)
Finally,

$$Q_{\mathrm{total}}=\frac{1}{1+{\displaystyle \frac{4\omega _0M^2}{\mu _0R}}}\left[\frac{1}{\pi }\left(\mathrm{log}\frac{4}{\gamma k_0r_0}1\right)+\frac{4M^2}{\mu _0R^2C}\right].$$
(44)
Naturally, one wants to minimize dissipation in the antenna in order to increase the antenna efficiency. However, from (44) it is seen that if the losses in the material become very small ($`R0`$), the total quality factor reduces to the well-known relation for an oscillatory circuit:

$$Q_{\mathrm{total}}=\frac{1}{\omega _0RC},$$
(45)
which means that most of the energy is stored in the medium layer and almost all the source power goes into heat (the antenna efficiency tends to zero).

Next, let us fix a given efficiency value fixing the resistance, the resonant frequency and the coupling. To decrease (44) one should increase the capacitance. In the limit $`C\mathrm{}`$. For a given fixed resonant frequency that means $`L0`$. One can see that the total impedance of a medium particle and the medium effective impedance become simply resistive in the limiting case:

$$Z_{\mathrm{medium}}=\frac{\omega _0^2M^2}{R}.$$
(46)
Hence, the deepest drop in the total quality factor can be achieved if the medium loading is equivalent to a simple resistor. The limiting value for the quality factor is for this case

$$Q_{\mathrm{total}}=\eta Q_{\mathrm{wire}},$$
(47)
as expected. Note also that when the inductance $`L`$ tends to zero, a better model should assume that also the mutual inductance $`M`$ tends to zero. In that case the effect due to the material cover simply disappears in the limit, as the particles forming the cover material are not excited.

Because the radiation resistance of the considered system does not depend on the layer parameters and, moreover, the covering does not affect the part of energy stored outside the cylinder, the antenna quality factor and the bandwidth associated with radiation does not change in the presence of covering. However, the total quality factor associated with the input impedance of the antenna as a circuit element is very much affected.

## 4 Discussion and conclusions

First of all, we see from Eq. (32) that the electric properties of the covering material have disappeared from the antenna quality factor expression. Such situation is usual when one considers thin coverings placed on top of electric conductors. But a more important conclusion following from Eq. (32) is that we cannot achieve any improvement in the antenna bandwidth even covering it with a (passive) magnetic material having a negative permeability. This follows from the causality limitations on the material parameters of any passive media with negligible losses :

$$\frac{(\omega ϵ_r)}{\omega }1,\frac{(\omega \mu _r)}{\omega }1.$$
(48)
Indeed, because of these limitations, the value of the first term in braces in Eq. (32) can only be increased compared to the case of empty filling (i.e. when $`\mu =\mu _0`$), in turn increasing $`Q`$ and narrowing the bandwidth.

Considering a dual problem involving a magnetic current line as the radiating source covered with the same material, one will get the same expression (32), only the magnetic permeability will be replaced by the electric permittivity. There will be no bandwidth improvement also.

In case of our example of an antenna with a dispersive and lossy covering, the additional stored reactive energy is of course also always positive, although the reactive part of the input impedance can cross zero and the derivative of $`X(\omega )`$ can change sign. The bandwidth can increase due to additional loss at a cost of reduced antenna efficiency. The formulas for the stored energy in this case are not universal and depend on the material structure, but the conclusion that the stored energy is positive in any passive system is universal.

To summarize, homogeneous coverings of radiating metal wires<sup>1</sup><sup>1</sup>1For finite-length wire antennas, the present theory is applicable to covers that do not substantially change the effective antenna height. by electrically thin passive material layers lead to reduced bandwidth or lower efficiency, whatever exotic properties these covering materials might have. The following alternative possibilities to improve antenna performance using complex materials can be identified:

* The use of radiating inclusions, and not really material coverings.
* Shells of resonant dimensions. New materials can offer more possibilities in optimizing resonant antennas. Apparently, this is the case considered in .
* Nonuniform coverings or material inclusions. This can modify the current distribution, possibly leading to increased bandwidth.
* Active materials. If the passivity requirement is dropped, it is in principle possible, for example, to realize wide-band weakly dispersive negative material parameters . The stored reactive energy can be negative (and bandwidth very large), meaning, effectively, that the whole volume filled by this material is the source of power.


text:
# Renormalization of the electron-phonon coupling in the one-band Hubbard model

## Abstract

We investigate the effect of electronic correlations on the coupling of electrons to Holstein phonons in the one-band Hubbard model. We calculate the static electron-phonon vertex within linear response of Kotliar-Ruckenstein slave-bosons in the paramagnetic saddle-point approximation. Within this approach the on-site Coulomb interaction $`U`$ strongly suppresses the coupling to Holstein phonons at low temperatures. Moreover the vertex function does not show particularly strong forward scattering. Going to larger temperatures $`kTt`$ we find that after an initial decrease with $`U`$, the electron-phonon coupling starts to increase with $`U`$, confirming a recent result of Cerruti, Cappelluti, and Pietronero. We show that this behavior is related to an unusual reentrant behavior from a phase separated to a paramagnetic state upon decreasing the temperature.

The relevance of phonons for high-temperature superconductivity has been debated since the discovery of the high-Tc cuprates. Recently, strong renormalization effects of the electrons near the Fermi surface, observed in angle-resolved photoemission in several cuprates, have been at least partially ascribed to phonons lanzara . Furthermore, quantum Monte Carlo simulations of the Hubbard-Holstein model suggest that the electron-phonon coupling shows forward scattering and no substantial suppression at large U and small dopings hanke , similar as in the $`1/N`$ expansion for the t-J model zeyher . On the other hand, it has been pointed out keller  that at small dopings the Kotliar-Ruckenstein (K-R) slave-boson approach KR  might yield results quite different from the 1/N expansion. Below we study the influence of strong electronic correlations on the electron-phonon coupling using the K-R approach. The quantity of interest is the static vertex function $`\mathrm{\Gamma }`$ which acts as a momentum-dependent, multiplicative renormalization factor for the bare electron-phonon coupling.

We consider the one-band Hubbard model on a square lattice with nearest and next-nearest neighbor hopping, $`t`$ and $`t^{}`$, respectively,

$$H=t\underset{i,j,\sigma }{}c_{j\sigma }^{}c_{i\sigma }^{}t^{}\underset{i,j,\sigma }{}c_{j\sigma }^{}c_{i\sigma }^{}+U\underset{i}{}n_in_i.$$
(1)
For the non-interacting system the dispersion relation is $`\epsilon _k=2t(\mathrm{cos}(k_x)+\mathrm{cos}(k_y))4t^{}\mathrm{cos}(k_x)\mathrm{cos}(k_y)`$, and the density of states has a logarithmic van Hove singularity at $`4t^{}`$. For this model we want to study the influence of electronic correlations on the coupling of electrons to an external field $`V_i`$. The bare coupling has the form $`H^{}=_{i,\sigma }n_{i\sigma }V_i`$. Writing $`V_i=gu_i`$, $`H^{}`$ describes also the interaction of electrons and atomic displacements $`u_i`$ with coupling constant $`g`$. The linear change in the one-particle Green’s function $`G(p)`$ due to $`V_q`$ is,

$$\frac{\delta G(p)}{\delta V_q}=G(p)\mathrm{\Gamma }(p,q)G(p+q),$$
(2)
with the charge or electron-phonon vertex $`\mathrm{\Gamma }(p,q)=\delta G^1(p)/\delta V_q`$. The components of the three-dimensional vectors $`p`$ and $`q`$ consist of a frequency and a two-dimensional momentum. For the calculation of the vertex we use the slave-boson technique of Kotliar and Ruckenstein KR . The basic idea of our approach sblin  is to calculate linear responses by linearizing the saddle-point equations for the perturbed system about the homogeneous saddle-point solution. We consider only paramagnetic solutions. Then there are three slave-bosons, $`e`$, $`p`$, and $`d`$, describing empty, singly, and doubly occupied sites, and two Lagrange parameters $`\lambda ^{(1)}`$ and $`\lambda ^{(2)}`$ enforcing consistency between slave-fermions and slave-bosons. The linear response to a charge-like perturbation of a given wave-vector can be determined by solving the $`5\times 5`$ system of linear equations given in sblin . Considering only static fields $`V_𝐪`$ we find

$$\mathrm{\Gamma }(𝐩,𝐪)=1+\frac{\delta \lambda ^{(2)}}{\delta V_𝐪}+z\left(\epsilon _𝐩+\epsilon _{𝐩+𝐪}\right)\frac{\delta z}{\delta V_𝐪}.$$
(3)
The first term in Eq. (3) is due to the explicit dependence of $`G^1`$ on $`V`$, the remaining terms are obtained by taking the derivative of the self-energy with respect to $`V`$. $`\mathrm{\Gamma }`$ does not depend on frequencies because we assumed zero frequency in $`q`$ and because the saddle-point self-energy is frequency-independent. $`z`$ is given by the Kotliar-Ruckenstein choice

$$z=\frac{(e+d)p}{\sqrt{1p^2d^2}\sqrt{1e^2p^2}}.$$
(4)
In the limit $`U\mathrm{}`$ our approach reduces to method (II) of Ref. keller , and we have checked that for large $`U`$ we recover the results given in their Fig. 1.

While it was shown that the slave-boson linear-response method gives very good results for the charge susceptibility (see, e.g., Fig. 1 of Ref. sblin  for a comparison with exact diagonalization), it is not clear a priori how well it will work for the charge vertex $`\mathrm{\Gamma }(𝐩,𝐪)`$. As a check we have calculated the static vertex for a small system using exact diagonalization. The result for the scattering of an electron from a state just below, to a state just above the Fermi surface is shown in Fig. 1. Considering the fact that in exact diagonalization the number of particles is fixed, while the slave-boson calculations are performed in the grand canonical ensemble, the agreement between both methods is remarkably good. This indicates that the slave-boson approach should work well at zero temperature.

To find out how well the slave-bosons work at finite temperatures we compare to the quantum Monte Carlo (QMC) calculations of Ref. hanke . Fig. 4 of that work shows the effective electron-phonon coupling $`g(𝐩,𝐪)`$, as defined in their Eq. (7), for the Hubbard model on an $`8\times 8`$ lattice with filling $`n=0.88`$, calculated at the lowest fermionic Matsubara frequency, for an inverse temperature of $`\beta =2`$. For comparison, we show in Fig. 2 the results of our slave-boson calculations for the same model at $`\omega =0`$ and a slightly different filling $`n=0.875`$. Also here we find a remarkable agreement. In particular, we find that after an initial decrease the coupling for forward scattering (small $`𝐪`$) starts to increase for $`U8`$. This seems to indicate that the slave-boson method also works well at finite temperatures. Moreover, Eq. (3) naturally explains why the QMC results for different electron momenta $`𝐩`$, shown in Fig. 4 of Ref. hanke , are so similar.

After comparing the results of the slave-boson linear response calculations to more accurate methods, which are, however, limited to small systems (exact diagonalization) or finite temperatures (quantum Monte Carlo), we now turn to very large systems at very low temperatures. First we calculate the electron-phonon vertex for electrons on the Fermi surface, where in Eq. (3) $`\epsilon _𝐩`$ and $`\epsilon _{𝐩+𝐪}`$ are both replaced by the Fermi energy of the non-interacting system. Fig. 3 shows the vertex for momentum transfer $`𝐪`$ along high symmetry lines in the Brillouin zone for the Hubbard model at essentially zero temperature. The effect of next-nearest neighbor hopping is illustrated by comparing calculations for $`t^{}=0`$ and $`t^{}=0.35t`$. We find that in both cases the on-site Coulomb interaction strongly reduces the electron-phonon coupling. This is not completely unexpected as the charge response should be strongly suppressed by an on-site Coulomb interaction. It is, however, in striking difference to the behavior at higher temperature (Fig. 2). Also, while $`\mathrm{\Gamma }(𝐪)`$ shows a broad peak around $`𝐪=0`$, we do not find a particularly pronounced forward scattering. In fact, the electron-phonon vertex is often strongest close to $`𝐪=(\pi ,0)`$. This is different from what was found within an $`1/N`$ expansion zeyher . The $`1/N`$ expansion relies on the smallness of $`1/\delta N`$, i.e., it breaks down at small dopings $`\delta `$. This can be seen from the fact that the charge-charge correlation function remains in leading order finite for $`\delta 0`$ though the exact correlation function vanishes in this limit. The Kotliar-Ruckenstein method, on the other hand, reproduces this limit correctly in leading order which makes it plausible that in this case the charge vertex is smaller than in the 1/N expansion, especially at smaller dopings keller . Which of the two methods is more reliable, in particular, near optimal doping, is not clear and can probably only be judged by comparison with exact numerical methods.

To assess the importance of the electron-phonon coupling for superconductivity we calculate the renormalization factor

$$\mathrm{\Lambda }_\alpha =\frac{_{\mathrm{FS}}\frac{dp}{|𝐯_𝐩|}_{\mathrm{FS}}\frac{dp^{}}{|𝐯_𝐩^{}|}g_\alpha (𝐩)\mathrm{\Gamma }(𝐩,𝐩^{}𝐩)g_\alpha (𝐩^{})}{z^2_{\mathrm{FS}}\frac{dp}{|𝐯_𝐩|}_{\mathrm{FS}}\frac{dp^{}}{|𝐯_𝐩^{}|}g_\alpha ^2(𝐩)}$$

for the pairing channels $`g_s(𝐩)=1`$, $`g_s^{}(𝐩)=\mathrm{cos}(p_x)+\mathrm{cos}(p_y)`$, $`g_{p_x}(𝐩)=\mathrm{sin}(p_x)`$, $`g_{d_{x^2y^2}}(𝐩)=\mathrm{cos}(p_x)\mathrm{cos}(p_y)`$, and $`g_{d_{xy}}(𝐩)=\mathrm{sin}(p_x)\mathrm{sin}(p_y)`$. $`\mathrm{\Lambda }_\alpha `$ is equal to the ratio $`\lambda _\alpha /\lambda _\alpha ^{(0)}`$, where $`\lambda _\alpha `$ and $`\lambda _\alpha ^{(0)}`$ denote the dimensionless electron-phonon coupling constants in the interacting and non-interacting cases, respectively. To judge the importance of forward scattering we also calculate the renormalization factor for transport,

$$\mathrm{\Lambda }_{\mathrm{tr}}=\frac{_{\mathrm{FS}}\frac{dp}{|𝐯_𝐩|}_{\mathrm{FS}}\frac{dp^{}}{|𝐯_𝐩^{}|}\mathrm{\Gamma }(𝐩,𝐩^{}𝐩)|𝐯(𝐩)𝐯(𝐩^{})|^2}{2z^2_{\mathrm{FS}}\frac{dp}{|𝐯_𝐩|}_{\mathrm{FS}}\frac{dp^{}}{|𝐯_𝐩^{}|}|𝐯(𝐩)|^2}.$$

The results are shown in Fig. 4. We find that for $`U10`$ the $`s`$-wave couplings decrease almost exponentially with $`U`$. For the special case of the Hubbard model with nearest neighbor hopping only ($`t^{}=0`$) we have $`\mathrm{\Lambda }_s^{}=\mathrm{\Lambda }_s`$, since $`g_s^{}`$ is constant on the Fermi surface. Moreover $`\mathrm{\Lambda }_{\mathrm{tr}}\mathrm{\Lambda }_s`$, reflecting that there is no pronounced forward scattering; only for larger $`U`$ does $`\mathrm{\Lambda }_{\mathrm{tr}}`$ become somewhat smaller than $`\mathrm{\Lambda }_s`$. But by then both coupling constants are already very small. The higher pairing channels are even weaker, starting from zero at $`U=0`$, going through a maximum around $`U=2`$ only to decay almost exponentially. We can thus conclude that within Kotliar-Ruckenstein slave-boson theory, restricting the system to be paramagnetic, the contribution of Holstein phonons to superconductivity should be very small.

We now come back to the surprising upturn of the electron-phonon vertex for $`U8`$ shown in Fig. 2 and also found in the QMC calculations of Ref. hanke . Calculating $`\mathrm{\Gamma }(𝐪)`$ at $`kTt`$, we indeed find a drastically different behavior than for $`T0`$: Instead of monotonically decreasing with $`U`$, the coupling starts to increase and develops a very strong forward scattering peak. An example is shown in Fig. 5. Looking at the charge response function $`\chi (𝐪)`$ in the paramagnetic phase we find that this behavior is a precursor of a phase-separation instability — a divergence of $`\chi (𝐪=0)`$. This has already been pointed out in Ref. emmanuele . Calculating the phase diagram, we find a very peculiar reentrant behavior around the phase separated region as shown in Fig. 6: When cooling down the system phase separates, but at low enough temperature it reverts back to the paramagnetic phase. Since in our calculations we only allow for a paramagnetic phase, other phases might mask the phase separation. Also, since slave-bosons may have problems at high temperatures KR , it is not clear if the Hubbard model really shows such an reentrant behavior. Nevertheless, we find a qualitatively similar behavior in the limit $`U\mathrm{}`$ in the gauge invariant $`1/N`$ expansion (i.e., a theory without Bose condensation). Phase separation at finite $`T`$ has also been proposed in Refs. Woelfle ; Imada ; Cappelluti . Moreover the good agreement with the quantum Monte Carlo calculations of Ref. hanke  suggests that our approach might indeed capture the relevant physics. It would therefore be interesting to test the phase diagram shown in Fig. 6 with QMC: A calculation for, e.g., $`\beta =1`$ and $`U=4\mathrm{}8`$, i.e., at temperatures and values of $`U`$, where QMC has little problems, should show clear signs of phase separation. Of course, these calculations should be done at $`\omega =0`$ as the extrapolation from finite Matsubara frequencies might be difficult close to an instability.

In conclusion, we have studied the influence of strong electronic correlations on the electron-phonon interaction for the Hubbard-Holstein model using the Kotliar-Ruckenstein slave boson method. For high temperatures the boundaries of the phase-separated state were determined in the $`TU`$ plane for different dopings and the increase of the static vertex $`\mathrm{\Gamma }`$ near the boundaries was studied, confirming and extending recent results of Ref. emmanuele . At low temperatures and moderate or small dopings we found that $`\mathrm{\Gamma }`$ does not exhibit pronounced forward scattering behavior and that $`\mathrm{\Gamma }`$ reduces dramatically the electron-phonon coupling. It seems that exact numerical calculations are necessary to judge the reliability of the $`1/N`$ and the Kotliar-Ruckenstein approaches.

We would like to thank O. Dolgov, O. Gunnarsson, W. Hanke, M. Lavagna, A. Muramatsu, and D. Vollhardt for fruitful discussions.


text:
# Rectangular-Mask Coronagraphs for High-Contrast Imaging

## 1 Introduction

Following our previous work on pupil-plane masks for high-contrast imaging in the context of terrestrial planet finding (see Kasdin et al. (2003); Vanderbei et al. (2003a, b)), we present in this paper further new pupil masks as well as some pupil-plane/image-plane combinations of masks. Most of our earlier pupil-plane mask designs have an inner working angle (IWA) of $`4\lambda /D`$, which implies a main mirror in the $`12`$ m range. One of the main objectives leading to the designs presented in this paper was to significantly decrease this inner working angle. If the inner working angle can be reduced to $`2\lambda /D`$, then the main mirror can be dropped down to the $`6`$ m range and still be capable of studying the same set of stars. This is our goal.

The simplest, but not best, pupil mask presented here achieves a contrast ratio of $`10^{10}`$ in four rectangular regions given by $`\{(\xi ,\zeta ):2<|\xi |<20,2<|\zeta |<20\}`$. The points of high-contrast that are closest to the center of the star’s image occur along the diagonals and therefore correspond to an inner working angle of $`2\sqrt{2}=2.8\lambda /D`$ (as usual, $`\lambda `$ is wavelength and $`D`$ is aperture).

The masks considered in this paper consist of rectangular arrays of rectangular opennings. They can be built as two identical striped masks one placed on top of the other but oriented so that the stripes on one mask are orthogonal to the stripes on the other (see Figure 1). We call the striped masks barcode masks. Each barcode mask need only provide a contrast ratio of $`10^5`$—it is then guaranteed that the two-mask overlay will achieve a contrast of $`10^{10}`$. This is important since it is quite feasible to check in a laboratory that a barcode mask achieves $`10^5`$ but is perhaps impossible to check that a mask achieves $`10^{10}`$ in any ground-based laboratory.

Furthermore, we show that one can get an even tighter IWA of $`1.4\sqrt{2}=2.0\lambda /D`$ by placing a rectangular pupil-plane mask in the Lyot plane of a traditional coronagraph and using an image-plane occulting mask consisting of a “plus-sign” shape, $`\{(\xi ,\zeta ):|\xi |<0.6,|\zeta |<0.6\}`$. We show that such a coronagraph can tolerate pointing errors as large as about $`0.05\lambda /D`$.

Finally, in the interest of producing an on-axis coronagraph, we present a Lyot-plane rectangular pupil mask which includes a $`2\%`$ 4-vane spider. A similar spider can then be placed in the entrance pupil and provide the opportunity to place a small secondary mirror at the center of the plane. Having such a spider over the entrance pupil provides the possibility of building an all-mirror on-axis system and in this way greatly reducing the negative effects of differential polarization that plague off-axis designs.

## 2 Barcode and Rectangular Pupil Masks

In this paper we assume that telescope optics follow the Fraunhofer approximation. Hence, given a pupil plane apodization function $`0A(x,y)1`$, the image-plane electric field is given by the two-dimensional Fourier transform

$$E(\xi ,\zeta )=\widehat{A}(\xi ,\zeta ):=e^{2\pi i(\xi x+\zeta y)}A(x,y)𝑑y𝑑x.$$
(1)
If the apodization function takes only the values zero and one, then the function represents a mask. If it depends only on one of the coordinates (either $`x`$ or $`y`$), then we say that it is a one-dimensional apodization or a one-dimensional mask. We also refer to one-dimensional masks as barcode masks. This family of masks were studied in Kasdin et al. (2004). In this paper, we are interested in masks that correspond to the tensor product of a pair of one-dimensional masks:

$$A=A_xA_yA(x,y)=A_x(x)A_y(y).$$
(2)
The electric field corresponding to a tensor product is itself a tensor product:

$$E=\widehat{A}=\widehat{A_xA_y}=\widehat{A_x}\widehat{A_y}.$$
(3)
In other words,

$$E(\xi ,\zeta )=\widehat{A_x}(\xi )\widehat{A_y}(\zeta ).$$
(4)
Tensor products of smooth apodizations were first proposed for terrestrial planet finding in Nisenson and Papaliolios (2001).

Since intensity is the square of the magnitude of the electric field, it follows that for a contrast of $`10^{10}`$, we seek apodization functions that provide the following contrast inequalities:

$$|E(\xi ,\zeta )|10^5E(0,0),(\xi ,\zeta )𝒪,$$
(5)
where $`𝒪`$ denotes the points in the image plane at which high contrast is to be achieved. Suppose as before that the apodization function is a tensor product. If the set $`𝒪`$ is a generalized rectangle

$$𝒪=\{(\xi ,\zeta ):\xi 𝒪_\xi ,\zeta 𝒪_\zeta \}$$
(6)
(i.e., $`𝒪=𝒪_\xi \times 𝒪_\zeta `$), then the contrast inequalities can be achieved by giving two one-dimensional apodizations that each achieve a contrast ratio of only $`10^5`$:

$`|\widehat{A_x}(\xi )|`$ $``$ $`10^{2.5}\widehat{A_x}(0),\xi 𝒪_\xi ,`$ (7)
$`|\widehat{A_y}(\zeta )|`$ $``$ $`10^{2.5}\widehat{A_y}(0),\zeta 𝒪_\zeta .`$ (8)
Figure 1 shows a rectangular mask corresponding to $`𝒪_\xi =𝒪_\zeta =\{\xi :2|\xi |25\}`$. The mask has a $`28.1\%`$ open area. Its Airy throughput, defined as

$$_{\xi _0}^{\xi _0}_{\xi _0}^{\xi _0}|E(\xi ,\zeta )|^2𝑑\zeta 𝑑\xi ,$$
(9)
where $`\xi _0`$ denotes the inner working angle for the barcode mask, is $`15.1\%`$.

Figure 2 shows another rectangular mask. For this mask, a $`2\%`$ central obstruction was imposed on the design. With such a central obstruction (and implied spiders) one can hang and hide a small secondary mirror and therefore build a telescope with an on-axis optical path. In this design, we have used $`𝒪_\xi =𝒪_\zeta =\{\xi :2|\xi |22\}`$. The mask has a $`4.6\%`$ open area. Its Airy throughput, is $`0.41\%`$, which is somewhat disappointing. In the Section 4, we present an on-axis design that achieves a much higher throughput.

We end this section by point out that PSFs associated with these masks depend on wavelength only in that the inner working angle is measured in units of $`\lambda /D`$, where $`\lambda `$ is wavelength and $`D`$ is aperture. Hence, measured in radians, at longer wavelengths the inner working angle is correspondingly larger.

## 3 A Lyot Coronagraph

Consider an imaging system consisting of an entrance pupil, a first image plane, a reimaged (Lyot) pupil plane, and a final image plane (containing an imaging device). We assume that each of the first three planes can be apodized/masked. Following Kuchner and Traub (2002), we let $`A`$ denote the mask function for the entrance pupil, $`\widehat{M}`$ the mask function for the first image plane, and $`L`$ the mask function for the Lyot pupil.

The electric field at point $`(\xi ,\zeta )`$ in the final image plane corresponding to an off-axis point source is a composition of mask multiplication and Fourier transform:

$$\left(L\left(\widehat{M}\left(AF_{\xi _0,\zeta _0}\right)\right)\right)(\xi ,\zeta )=\widehat{L}\left(\widehat{M}(\widehat{A}\delta _{\xi _0,\zeta _0})\right)(\xi ,\zeta ),$$
(10)
where $`F_{\xi _0,\zeta _0}`$ denotes the electric field at the entrance pupil corresponding to a point source from direction $`(\xi _0,\zeta _0)`$, $`\delta _{\xi _0,\zeta _0}`$ denotes a unit mass delta function at $`(\xi _0,\zeta _0)`$, hats denote Fourier transforms and stars denote convolutions.

It is easy to see that if each of $`L`$, $`M`$, and $`A`$ are tensor products ($`L=L_xL_y`$, $`M=M_xM_y`$, $`A=A_xA_y`$), then the electric field again factors into a product of electric fields:

$$\widehat{L}\left(\widehat{M}(\widehat{A}\delta _{\xi _0,\zeta _0})\right)(\xi ,\zeta )=\widehat{L_x}\left(\widehat{M_x}(\widehat{A_x}\delta _{\xi _0})\right)(\xi )\widehat{L_y}\left(\widehat{M_y}(\widehat{A_y}\delta _{\zeta _0})\right)(\zeta )$$
(11)
We use this factorization to greatly simplify the coronagraph optimization problems described in the next section.

## 4 Lyot Coronagraph Optimization Problem and Numerical Results

As in our previous papers (Kasdin et al. (2003); Vanderbei et al. (2003a, b)), we use numerical optimization to maximize some measure $`\theta `$ of throughput subject to imposed contrast constraints. For the designs presented in this section, we solve the following one-dimensional barcode optimization problem:

$$\begin{array}{ccccc}\hfill \text{maximize }& \hfill \theta & & & \\ \hfill \text{subject to }& \hfill \theta & \widehat{L}\left(\widehat{M}(\widehat{A}\delta _{\xi ,\zeta })\right)(\xi ,\zeta ),\hfill & (\xi ,\zeta )𝒪,\hfill & \\ & \hfill \left|\widehat{L}\left(\widehat{M}(\widehat{A}\delta _0)\right)(\xi ,\zeta )\right|& \frac{\widehat{L}\left(\widehat{M}(\widehat{A}\delta _{\xi ,\zeta })\right)(\xi ,\zeta )}{10^{2.5}},\hfill & (\xi ,\zeta )𝒪.\hfill & \end{array}$$
(12)
The first set of constraints say that $`\theta `$ is a lower bound on the electric field of a unit off-axis source (i.e, a planet as bright as a star) as the off-axis source varies over the high-contrast region. The second set of constraints say that the magnitude of the electric field at $`\xi 𝒪`$ due to an on-axis star is smaller than the electric field at the same point in the image plane due to a much fainter planet whose image is centered at this point $`\xi `$. The functions $`L`$, $`M`$, and $`A`$ are functions of a single real variable and represent either the $`x`$ or the $`y`$ component of a tensor product that one forms to make the final two dimensional image. By specifying prescribed functions for $`A`$ and $`M`$ (actually, $`\widehat{M}`$) we can optimize only the function $`L`$. In this case, the problem is an infinite dimensional linear programming problem. Discretizing the image and pupil planes reduces the problem to a finite dimensional linear programming problem that can be solved numerically.

We now give two specific examples of Lyot coronagraphs. For the first example, we specify that

$`A(x)=\{\begin{array}{cc}1\hfill & |x|\frac{1}{2}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (15)
$`\widehat{M}(\xi )=\{\begin{array}{cc}1\hfill & 0.6|\xi |84\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (18)
$`|L(x)|\{\begin{array}{cc}1\hfill & |x|\frac{1}{2}\hfill \\ =0\hfill & \text{otherwise}\hfill \end{array}`$ (21)
$`𝒪=\{\xi :1.4|\xi |21\}`$ (22)
The optimal barcode mask is shown in Figure 3 together with the corresponding rectangular mask.

The second example is the same as the first except that we impose a central obstruction in the pupil-plane masks:

$`A(x)=\{\begin{array}{cc}1\hfill & 0.01|x|\frac{1}{2}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}`$ (25)
$`|L(x)|\{\begin{array}{cc}1\hfill & 0.01|x|\frac{1}{2}\hfill \\ =0\hfill & \text{otherwise}\hfill \end{array}`$ (28)
The optimal barcode mask is shown in Figure 4 together with the corresponding rectangular mask.

Since the designs presented in this section involve an image plane mask, the associated PSFs depend on wavelength in a more complicated way than the designs given in Section 2. In fact, as Figure 5 shows, the loss of contrast occurs at wavelengths shorter than about $`90\%`$ and longer than $`105\%`$ of the design point.

## 5 Conclusions

In this paper we have presented some new masks for high-contrast imaging. The best design

1. achieves an inner working angle of $`2.0\lambda /D`$ along the four diagonal directions,
2. includes a four-vane spider so that it can be used in an on-axis telescope design, and
3. consists of two orthogonally oriented barcode masks each one, by itself, being designed to provide only $`10^5`$ level of contrast.

We include an analysis of the effects of pointing error and chromaticity.

Acknowledgements. We are grateful for the continuing support we have received for this work from NASA and JPL. The first author also wishes to acknowledge support from the NSF (CCR-0098040) and the ONR (N00014-98-1-0036).


text:
# Structural Vulnerability of the North American Power Grid

## Abstract

The magnitude of the August 2003 blackout affecting the United States has put the challenges of energy transmission and distribution into limelight. Despite all the interest and concerted effort, the complexity and interconnectivity of the electric infrastructure have so far precluded us from understanding why certain events happened. In this paper we study the power grid from a network perspective and determine its ability to transfer power between generators and consumers when certain nodes are disrupted. We find that the power grid is robust to most perturbations, yet disturbances affecting key transmision substations greatly reduce its ability to function. We emphasize that the global properties of the underlying network must be understood as they greatly affect local behavior.

During the past decades the North American power infrastructure has evolved into what many experts consider the largest and most complex system of the technological age. Geographically, the power grid forms a network of over 1 million kilometers of high voltage lines that are continuously regulated by sophisticated flow control equipmentroadmap . As a result of the recent deregulation of power generation and transmission, about one-half of all domestic generation is now sold over ever-increasing distances on the wholesale market before it is delivered to customersroadmap . Consequently the power grid is witnessing power flows in unprecedented magnitudes and directionsrely .

As the power grid increases in size and complexity, it is becoming more important to understand the emergent behaviors that can take place in the system. Performing an analytic description of the electromagnetic processes integrated over the whole grid is a daunting, if not impossible, task. Instead the power industry must resort to constructing models that can be used to simulate the network’s response to various external parameters. Generally these models attempt to simulate actual electrical flow characteristics in smaller systems like a single distribution gridflow . In the present analysis we propose an alternative approach based on recent advances in understanding the structure of large complex networksab02 . We choose to investigate the network representation of the power grid from a topological perspective with the hope of finding properties and behaviors that transcend the abstraction.

We have built the network model based on data stored in the POWERmap mapping system developed by Plattsplatts , the energy information and market services unit of the McGraw-Hill Companies. This mapping system contains information about every power plant, major substation and $`115765`$ kV power line of the North American power grid. Our model represents the power grid as a network of 14,099 nodes (substations) and 19,657 edges (transmission lines). We distinguish three types of substations: generators are the sources for power, transmission substations transfer the power among high voltage transmission lines, and distribution substations are at the “outer edge” of the transmission grid, and the centers of local distribution grids. Only the identity of generating substations was directly available from our data sources. We identify distribution substations by the criterion of having a single high-voltage transmission line connected to them, with the expectation that the flow out of them is continued on smaller voltage feeder lines leading to consumersmiss . A total of 1633 nodes are power plants, we classify 2179 nodes as distributing substations, with the rest being labeled as transmission substations.

We consider the power from a generator to be accessible to a consumer if there is a path of transmission lines between the two. In practice, the existence of a connection between two substations does not always imply that power can be transferred across it as there may be capacity or other constrains present. By ignoring these our model provides an idealized view, a “best case scenario” regarding the characteristics of the grid. We find that the network representation of the power grid contains a single connected component, meaning that there is a path of transmission lines between any power plant and any distribution substation. This observation implies that in the best case scenario each distribution substation can possibly receive power from any generator.

Recent advances in mapping the topology of complex networks have uncovered that a large fraction of them are highly heterogeneous with respect to the number of edges incident on a node (also called the node degree). In these networks the majority of the nodes have low degrees, but there is a continuous hierarchy of high-degree nodes (hubs) that play an important role in the system. The degree distribution of these networks follows a power-law $`P(k)k^\gamma `$ with the exponent $`\gamma `$ mostly between $`2`$ and $`3`$. It was demonstrated both numerically and analytically that these so-called scale-free networks are resilient to the random loss of nodes, but are vulnerable to attacks targeting the high-degree hubsajb00 ; cnsw00 ; cebh01 . Therefore it is important both from a theoretical and practical standpoint to determine whether the connectivity of the power grid is reliant on a small set of hubs and whether their loss will cause a large-scale breakdown of the power grid’s transmission capability.

As the node degree is a good indicator of its topological importance, we first determine the degree distribution of the power grid. We find that the cumulative degree distribution defined as $`P(k>K)=_{k>K}P(k)`$ follows an exponential

$$P(k>K)exp(0.5K)$$
(1)
(see Fig. 1). This functional form agrees with previous results on the degree distribution of the Western power gridasbs00  and its classification as a single-scale network. The cumulative degree distribution shows that the probability of high-degree nodes is less than in a scale-free network, but higher than in a random network with the same number of nodes and edges. Power engineering principles suggest that the hubs of the power grid should belong to central station generators, and transmission substations should not have more than a few edges. Indeed, the inset to Fig. 1 shows that the fraction of generating substations among substations of a given degree increases with this degree. Surprisingly, however, there are several high-degree transmission substations (e.g. $`50`$ have degree higher than $`10`$), including the node with highest degree.

As the role of the power grid is to transport power from generators to consumers, a possible measure for the importance of a node corresponding to a substation is its betweenness (or load)gkk01 ; n01 . The betweenness of a node in a network is defined as the number of shortest paths that traverse itgkk01 ; n01 . Assuming that power is routed through the most direct path, the betweenness of a substation is a proxy for how much power it is transmitting, and for this reason we will use the alternative term load to denote it. Since it is the transmission substations’ role to route power from generators to distribution substations, we focus our attention to them. We determine the shortest paths starting from all generation substations and ending on an all other reachable substations. For each transmission node we accumulate the number of paths that pass through it; being at the start or at the end of a path does not count. The highest possible load is $`1633\times 1246620`$ million. We find that substations can have a load anywhere between 1 and 4 million, and determine the cumulative load distribution, i.e. the probability that a node’s load $`l`$ is larger than a given value $`L`$ (see Fig. 2). The functional form of the cumulative load distribution is

$$P(l>L)(2500+L)^{0.7}$$
(2)
Fig. 2 illustrates that $`40\%`$ of the substations participate in tens or hundreds of paths only, but $`1\%`$ of them are part of a million or more paths. These high-load substations, although possibly not hubs regarding their degree, play an important role in power transmission.

A fundamental requirement of the power grid is robustness, the ability to withstand and tolerate errors (random failure) and targeted attacksajb00 ; cebh01 ; cnsw00 . To ensure the reliability of power distribution, the transmission grid was conceived in such a way that there is more than one electrical path between any two points in the system power\_book . We wanted to verify whether the actual topology of the current power grid has this feature of global redundancy, or it has lost it during its growth and evolution. A possible measure of network redundancy is the so-called edge range, defined as the distance between the two endpoints of an edge if the edge connecting them were removed mnl02 . The inset of Fig. 2 shows the frequency of different edge ranges $`r`$ plotted as a function of $`r^1`$. We find that parallel edges and short alternative paths are fairly frequent. However, around $`15\%`$ of the edges in the power grid have an infinite range. In addition to the $`2179`$ edges ending in distribution substations, close to $`900`$ edges connecting generators and/or transmission substations are radial. These radial edges represent a clear vulnerability, as their loss disconnects their endpoints and creates isolated clusters in the power grid.

While the connectedness of the power grid allows for the transmission of power over large distances, it also implies that local disturbances propagate over the whole grid. The failure of a power line due to lightning strike or short-circuit leads to the overloading of parallel and nearby lines. Power lines are guarded by automatic devices that take them out of service when the voltage on them is too high. Generating substations are designed to switch off if their power cannot be transmitted; this protective measure has the unwanted effect of diminishing power for all consumers. Another possible consequence of power line failure is the incapacitation of transmission substations, possibly causing that the power from generators cannot reach distribution substations and ultimately consumers.

In the unperturbed state each distribution substation can receive power from any of the $`N_g=1633`$ generators. As substations lose function, the number of generators connected to (and able to feed) a certain distribution substation $`i`$, $`N_g^i`$, decreases. We introduce the concept of connectivity loss to quantify the average decrease in the number of generators connected to a distributing substation,

$$CL=1\frac{N_g^i}{N_g}_i,$$
(3)
where the averaging is done over every distributing substation. In summary, the connectivity loss measures the decrease of the ability of distribution substations to receive power from the generators, and in the following we will express it as a percentage.

First we investigate the effect that the failure of a power-generating substation has on consumers. Since initially the network contains a single connected component every consumer can reach all generators, and their connectivity is $`100\%`$. As the number of generators decreases this value will decrease due to both loss of the generators themselves and due to loss of routing capabilities at the generating substation level. We remove nodes corresponding to generators either randomly, or in the decreasing order of their degrees, and monitor the connectivity loss as a function of the fraction of generators missing. The minimum possible loss is equal to the fraction $`f_g`$ of inactive generators and is due to the loss in generation only (straight line on Fig. 3). We find that the connectivity loss caused by removing power substations remains very close to this minimum value (Fig. 3), even though generating substations tend to be the largest hubs in the system. The removal of generating substations does not alter the overall connectivity of the grid thanks to a high level of redundancy at the power generating substation level.

The situation can be dramatically different when the nodes that we remove are transmission nodes. If the power grid were highly redundant the loss of a small number of transmission substations should not cause power loss as power is rerouted through alternative paths. We find that even the removal of a single transmission node causes a slight connectivity loss. We remove transmission nodes one by one, first randomly, then in the decreasing order of their degree or load. For a random failure the connectivity loss is fairly low and stays proportional with the number of nodes lost. The connectivity loss is significantly higher, however, when targeting high degree or high load transmission hubs (Fig. 4). The grid can withstand only a few failures of this nature before considerable parts of the network become disconnected leading to substantial connectivity loss at consumer level. For example, failure of only $`4\%`$ of the nodes with high load may cause up to $`60\%`$ loss of connectivity. We also study an algorithm where we periodically recalculate the load of all transmission nodes during node removal, and select the nodes with highest load to be deleted next. This is a possible illustration of a propagating (cascading) power failure, where it is more likely that substations that have the highest load in the perturbed configuration will fail next. Fig. 4 illustrates that this cascading failure has the most damaging effect, as the loss of only $`2\%`$ of the high-load transmission substations leads to a connectivity loss of almost $`60\%`$, and all distribution substations become virtually powerless at $`f_t8\%`$. In conclusion, the transmission hubs ensuring the connectivity of the power grid are also its largest liability in case of power breakdowns.

This vulnerability of the electric power grid is inherent to its organization and therefore cannot be easily addressed without significant investment. Possible solutions include increasing the redundancy and capacity of the existent structure or decreasing the reliance on transmission by incorporating more generation at the distribution substation level. Such distributed generation by small local plants can supplement power from the grid under normal operation conditions and can greatly mitigate the effects of blackouts on the population. Targeted use of generation located near the point of use might prove to be the only viable economical alternative.

###### Acknowledgements.

The authors wish to thank Donna Heimiller and Steven Englebretson for their help in obtaining the POWERmap network data. This research was partially supported by the Midwest Research Institute (contract number AAX-3–33641-01).


text:
# Waveguiding properties of optical vortex solitons

## I Introduction

Vortices are fundamental localized objects which appear in many branches of physics, ranging from the physics of liquid crystals to the dynamics of superfluids and Bose-Einstein condensates pismen . Different types of vortices can be found and identified in optics; one of the simplest objects of this kind is a phase singularity in an optical wave front which is associated with a phase dislocation carried by a diffracting optical beam review .

When such an optical vortex propagates in self-defocusing nonlinear medium, the vortex core with a phase dislocation becomes self-trapped, and the resulting stationary singular beam is known as an optical vortex soliton swartz . Such optical vortex solitons exist on a nonvanishing background wave, and they represent a two-dimensional generalization of the so-called dark solitons book . They can be generated in experiment as self-trapped singularities of broad beams by using a holographic phase mask, as was demonstrated for different types of nonlinear defocusing media exp1 ; exp2 ; exp3 . As a matter of fact, this type of optical vortex solitons demonstrates many common features with the vortices observed in superfluids and Bose-Einstein condensates bec .

In self-focusing nonlinear media, optical vortices can exist as ring-like optical beams with zero intensity at the center carrying a phase singularity kruglov . However, due to the self-focusing nature of nonlinearity such ring-like vortex beams become unstable, and they are known to decay into several fundamental optical solitons flying off the main ring firth . This effect was observed experimentally in different nonlinear systems, including the saturable Kerr-like nonlinear media, biased photorefractive crystals, and quadratic nonlinear media in the self-focusing regime (see details and references in Ref. book ).

Waveguides induced by optical vortices in both linear and nonlinear regimes are of a special interest because this type of waveguides is robust and can be made reconfigurable. Earlier numerical and experimental results wave1 ; wave2 ; wave3  indicate that many of the vortex waveguiding properties should be similar to those of planar dark solitons book . Moreover, the vortex-induced waveguides can guide large-amplitude beams beyond the applicability limits of the linear guided-wave theory, and, together with the guided beam, they can form a special type of the vortex-mode vector soliton or its various generalizations vortex ; dipole ; dipole2 ; dipole3 .

In this paper, we carry out a systematic analysis of the waveguiding properties of the vortex solitons and vortex-mode vector solitons in saturable nonlinear media, for both self-defocusing and self-focusing nonlinear media. We follow the earlier analysis of Ref. wave3  and examine two major regimes of the vortex waveguiding. In the linear regime, the guided beam is weak, and the analysis of the vortex waveguiding properties is possible by approximate analytical methods, reducing the problem to a well–known analysis of the linear guided-wave theory (as an example, see Ref. wave2 ). The most interesting nonlinear regime corresponds to large intensities of the guided beam, and it gives rise to composite (or vector) solitons with a vortex component, that should be identified and analyzed numerically.

We also analyze the vortex waveguiding and vortex-mode vector solitons in a self-focusing nonlinear medium where, as is well known, the vortex beam becomes self-trapped and it also displays strong azimuthal modulational instability. However, we demonstrate that a mutual incoherent coupling between the vortex waveguide and a large-amplitude guided mode it guides can provide a strong stabilizing mechanism for stable two-component vortex solitons to exist in such media.

The paper is organized as follows. The next section (Sec. II) presents our model which takes into account both the incoherent coupling between the waveguide and the localized mode it guides, and the saturable nature of nonlinearity usually realized in experiment. The model is valid for both self-defocusing and self-focusing media. Section III is devoted to the study of the properties of optical vortex solitons and their guiding properties in self-defocusing media, whereas Sec. IV summarizes our results for the vortex waveguides and vector solitons in self-focusing nonlinear media. In particular, for the first time to our knowledge, we demonstrate that the vortex-mode vector solitons formed by the vortex waveguide together with a large-amplitude mode it guides can show the bistability property, when two types of different composite soliton structures exist for the same value of the soliton power. We also discuss, by means of direct numerical simulations, the stability properties of the vortex-mode vector solitons in a self-focusing nonlinear medium. In particular, we reveal an effective method of suppression of the azimuthal instability through the stabilization effect that the guided mode exerts on the hosting vortex, the latter is known to be unstable by itself in focusing media. Finally, Sec. V concludes our paper.

## II Model

In order to study the guided modes of a vortex-induced waveguide in both self-focusing and self-defocusing media, we consider the interaction of two mutually incoherent optical beams propagating in a nonlinear saturable medium. The equations of motion for two beams can be presented in the dimensionless form as follows,

$$\begin{array}{c}i\frac{u}{z}+\mathrm{\Delta }_{}u+\eta \frac{(|u|^2+\mu |v|^2)u}{1+\sigma (|u|^2+|v|^2)}=0,\hfill \\ i\frac{v}{z}+\mathrm{\Delta }_{}v+\eta \frac{(|v|^2+\mu |u|^2)v}{1+\sigma (|u|^2+|v|^2)}=0,\hfill \end{array}$$
(1)
where $`u`$ and $`v`$ are the dimensionless amplitudes of two fields, the parameter $`\sigma `$ characterizes the nonlinearity saturation effect, and the incoherent mode interaction is described by the coupling parameter $`\mu `$. The spatial coordinate $`z`$ is the propagation direction of the beams, and $`\mathrm{\Delta }_{}`$ stands for the transversal part of the Laplacian operator in the cylindrical coordinates $`r=(x^2+y^2)^{1/2}`$ and $`\varphi =\mathrm{tan}^1(y/x)`$. The sign parameter $`\eta =\pm 1`$ defines the type of the nonlinear medium under consideration, which is self-focusing, for $`\eta =+1`$, or self-defocusing, for $`\eta =1`$, respectively.

The model (1) provides a straightforward generalization to a number of important cases studied earlier. In particular, the limit $`\sigma 0`$ corresponds to the Kerr medium discussed, for example, in Refs. wave1 ; wave3 ; haelt , whereas the saturable case at $`\mu =1`$ corresponds to the incoherent beam interaction in photorefractive nonlinear media where different types of composite vector solitons have been predicted theoretically and observed in experiment vortex ; dipole ; dipole2 ; dipole3 . In particular, some examples of the two-component vector solitons composed of a vortex component and the localized modes it guides were presented earlier by Haelterman and Sheppard haelt  for the defocusing Kerr nonlinearity.

We look for stationary solutions of the system 1 in the form of a radially symmetric single-charged vortex beam in the field $`u`$ of the form

$$u(r,\varphi ;z)=u(r)e^{i\varphi }e^{i\gamma z},$$
(2)
where $`\gamma `$ is the propagation constant of the vortex mode, which if considered normalized takes the values +1 or -1 in the focusing and defocusing cases respectively, and consequently we have $`\gamma =\eta `$. For $`r\mathrm{}`$, the amplitude $`u(r)`$ approaches a constant value, for the self-defocusing medium, or zero, for the self-focusing case, respectively. We assume that the vortex guides the second beam given by the expression

$$v(r,\varphi ;z)=v(r)e^{il\varphi }e^{i\beta z},$$
(3)
where $`\beta `$ is the propagation constant of the guided beam. In the expressions (2), (3) the functions $`u(r)`$ and $`v(r)`$ are the radial envelopes of the interacting fields, and $`l`$ is the angular momentum of the guided mode.

In this way, the following set of stationary $`z`$-independent equations are obtained:

$$\begin{array}{c}\eta u+\mathrm{\Delta }_ru\frac{1}{r^2}u+\eta \frac{(u^2+\mu v^2)u}{1+\sigma (u^2+v^2)}=0,\hfill \\ \beta v+\mathrm{\Delta }_rv\frac{l^2}{r^2}v+\eta \frac{(v^2+\mu u^2)v}{1+\sigma (u^2+v^2)}=0,\hfill \end{array}$$
(4)
where $`\mathrm{\Delta }_r(1/r)d/dr(rd/dr)`$. The model 4 describes different types of localized solutions for the vortex and the localized modes it guides, and it has been further analyzed by the numerical relaxation methods to classify its localized solutions for both defocusing and focusing cases, respectively.

In the defocusing case when $`\eta =1`$, we have $`uu_0`$ and $`v0`$ for $`r\mathrm{}`$. Using this asymptotic conditions, we obtain

$$\left(1\frac{u_0^2}{1+\sigma u_0^2}\right)u_0=0,$$
(5)
and find the result for the vortex background amplitude $`u_0`$ in the form, $`u_0=(1\sigma )^{1/2}`$.

Using the second equation for $`v`$, and neglecting the terms which most quickly tend to zero when $`r\mathrm{}`$, we obtain:

$$\beta v+\frac{d^2v}{dr^2}\frac{\mu u_0^2}{1+\sigma u_0^2}v=0.$$
(6)
In order to have the exponentially decaying function $`v`$ (for $`r\mathrm{}`$), the following condition should be satisfied,

$$\beta +\frac{\mu u_0^2}{1+\sigma u_0^2}<0.$$
(7)
that gives as the result

$$\beta +\mu >0\beta >\mu ,$$
(8)
which defines the existence domain threshold as the line $`\beta =\mu `$. This result is valid for both fundamental and first-order guided modes.

## III Vortex guided modes

First, we consider the case of self-defocusing nonlinearity when $`\eta =1`$. In this case, the vortex amplitude $`u(r)`$ approaches a finite-amplitude asymptotic value $`u_0`$ for $`r\mathrm{}`$, where $`u_0`$ can be found from the asymptotic analysis, $`u_0=(1\sigma )^{1/2}`$. With the help of numerical relaxation methods, we find various types of localized modes that describe a vortex and the localized guided mode with the amplitude $`v(r)`$ satisfying the condition $`v(r)0`$ for $`r\mathrm{}`$. This condition requires that $`\beta <0`$, since the field $`v(r)`$ is localized due to the presence of the vortex waveguide only, and it describes a guided mode.

Figure 1 shows several characteristic solutions of this type for particular values of $`\mu `$ and $`\sigma =0.5`$, whereas the existence domain for the solutions with the fundamental guided mode (i.e. when $`l=0`$) is shown in Fig. 2. For the values of the propagation constant $`\beta `$ close to the cutoff ($`\beta _c0.79`$ at $`\mu =1`$), the second field presents a low-amplitude guided mode that can be analyzed by means of the linear guided wave theory. However, for values of $`\beta `$ close to the threshold, the amplitude of the guided mode becomes comparable to that of the vortex and its width grows, as shown for the case E in Fig. 1. In this latter case, the vortex is deformed as well, so the resulting structure describes a kind of a two-component vector soliton with the vortex-like mode.

In Fig. 3, we show the existence domain for the first vortex-like ($`l=1`$) mode of the induced vortex waveguide. Similar to the case of the fundamental guided mode, the linear theory is applicable near the mode cut-off, as is shown for the mode C, where the mode grows and it deformes the vortex waveguide substantially, see the case A near the upper boundary of the existence domain in Fig. 3 (top).

## IV Self-trapped vortex beams

### IV.1 Vortex-mode vector solitons

In the case of the self-focusing nonlinear medium, we put $`\eta =+1`$, and consider the self-trapped vortex beams with the asymptotic behavior $`u(r)0`$ as $`r\mathrm{}`$. Besides that, it is required that $`v(r)0`$ and $`\beta >0`$, so that $`v(r)`$ describes a localized mode. At the origin ($`r=0`$), the boundary condition for the vortex is $`u=0`$ and, if we seek the solutions corresponding to the fundamental guided mode ($`l=0`$), we have $`v=v_0`$ (or $`dv/dr=0`$) for the field $`v(r)`$ at $`r=0`$.

In Fig. 4, we show some examples of the two-component localized solutions which describes a fundamental (no nodes) beam guided by the self-trapped vortex in the form of a vortex-mode soliton. The existence domain for this kind of solutions is shown in Fig. 5 on the parameter plane ($`\mu `$, $`\beta `$), calculated numerically for $`\sigma =0.5`$. Close to the cutoff, the guided mode has a small amplitude (see the cases A and D2), and the vortex is only weakly distorted. However, when the propagation constant $`\beta `$ is close to the upper threshold, the field $`v(r)`$ grows in the amplitude affecting strongly the vortex mode.

In order to describe the bistable vector solitons, in Fig. 6 we plot the bifurcation diagram of the two-component localized solutions, for the partial and total beam powers. The power of the composite vortex-mode solitons $`P_{\mathrm{total}}`$ originates at the bifurcation point $`O_1`$ where the mode $`v`$ is small and it can be described by the linear theory. For larger value of $`\beta `$, the curve bends, and then it merges with the other partial power curve $`P_v`$ at the bifurcation point $`O_2`$ (see Fig. 6). The bistable solutions B1 and B2 correspond to a single value of the propagation constant $`\beta `$ in the bistability domain. Importantly, two solutions have different stability properties, and only one of the solutions is stable, as shown in Fig. 7.

### IV.2 Vortex stabilization

The incoherent interaction between the vortex and the localized mode it guides has the character of attraction, and it may provide an effective physical mechanism for stabilizing the vortex beam in a self-focusing nonlinear medium. Indeed, it is well-known that the vortex beam becomes unstable in a self-focusing nonlinear medium due to the effect of the azimuthal modulational instability. In this case, the vortex splits into the fundamental beams that fly off the main vortex ring firth . On the other hand, the bright solitons are known to be stable in such media. We expect that a mutual attraction of the components in a two-component system would lead to a counter-balance of the vortex instability by the bright component if the amplitude of the latter is large enough.

To carry out this study, we consider a two-component composite structure consisting of a vortex beam together with the fundamental mode of the waveguide it guides, described by Eqs. 1 (at $`\eta =1`$ and $`\mu =1.0`$). To study the mode stability, we propagate the stationary soliton solutions numerically.

In Fig. 8, we compare the vortex stability for the scalar and two-component systems. In the top row, we show the propagation of a vortex alone in the scalar model; the vortex breaks into two solitons which fly away after some propagation distance. In the bottom row, we show the propagation of two coupled components (the vortex and bright mode it guides). Due to a strong coupling between the modes, the propagation of the vortex is stabilized and its decay can be delayed dramatically, as shown in Fig. 8(lower row), or become completely stable, similar to the case shown in Fig. 7(low row). We confirm this stabilization mechanism by performing a similar study for a Gaussian input beam of the bright component instead of the exact stationary state, as is easier to realize in experiment.

## V Conclusions

We have studied the properties of the vortex-induced optical waveguides in both self-focusing and self-defocusing saturable nonlinear media. We have calculated numerically the existence domains of the composite solitons created by an optical vortex and the localized mode it guides, and have analyzed the properties of the vortex waveguides in both linear and nonlinear regimes. In the case of self-defocusing nonlinear media, we have identified the regimes when the vortex soliton can guide both the fundamental and first-order guided modes creating an effective multi-mode waveguide. In the case of self-focusing nonlinear media, we have studied the effect of the strong coupling between the vortex waveguide and the guided mode on the vortex stability. In particular, we have described a novel mechanism of stabilizing the vortex azimuthal modulational instability by a large-amplitude guide mode, and we have revealed the existence of optical bistability for the vortex-mode composite solitons. We believe similar results can be obtained for vortices in other types of nonlinear models, and they can be useful for other fields such as the physics of the Bose-Einstein condensates of ultra-cold atoms.

## Acknowledgements

We thank Sergey Odoulov and Mikhail Vasnetsov for the invitation to submit our contribution to the special issue, and we dedicate this article to Professor Marat Soskin, a pioneer of the physics of optical vortices and singular optics, on the occasion of his 75-th birthday.

This work was partially supported by the Australian Research Council. J.R. Salgueiro acknowledges a posdoctoral fellowship granted by the Secretaría de Estado de Educación y Universidades of Spain and partially supported by the European Social Fund.