text:
# Thin Film Diblock Copolymers in Electric Field: Transition from Perpendicular to Parallel Lamellae

## I Introduction

Diblock copolymers are known to self-assemble into a variety of ordered structures, with a length scale ranging from nanometers to micrometers These phases have potential applications in nanolithographic templates,$`^{\text{[1]}}`$ waveguides $`^{\text{[2]}}`$ and dielectric mirrors.$`^{\text{[3]}}`$ The length scale and morphology in the melt can easily be adjusted by controlling the fraction $`f=N_A/N`$ of A monomers in an A/B chain of $`N=N_A+N_B`$ monomers, and the temperature $`T`$.$`^{\text{[4, 5, 6]}}`$

In experiments one often encounters samples in the lamellar phase (made up of alternating planar A- and B-rich domains) in which the melt is only partially ordered. Ordered microdomains, or grains (typical size in the micrometer range), with grain boundaries between them are defects that cost energy. In order to anneal these defects, and to create a perfect alignment of the lamellae, several techniques have been used. In the bulk, mechanical shear proved to be a successful technique. Alignment by application of an electric field $`^{\text{[7, 8, 9, 10]}}`$ is also possible, but for macroscopic samples it requires high voltage difference between the two bounding electrodes. Nevertheless, this technique is especially suitable for thin films, because the thickness involved makes the required large fields (typically $`1030\mathrm{V}/\mu \mathrm{m}`$) accessible.

We consider in this paper thin films of lamellar diblock copolymers, under the influence of a perpendicular electric field. Initially, the lamellae are parallel to the confining surfaces, because of preferential short-range interactions with the surfaces. In section II we consider diblock copolymers in the weak segregation regime. We show in section III that electric field applied perpendicular to the surfaces can cause the melt to transform from a parallel to a perpendicular orientation through a first-order phase transition. The critical field $`E_c`$ for this transition is caused by a competition between the electric field and surface interactions. In section IV we investigate the thin-film alignment for diblocks in the strong segregation regime. In this regime, the surface correlations are finite, and thus the range of parallel ordering induced by the surfaces is finite as well. We give the transitions between parallel, perpendicular and mixed lamellae in terms of the system parameters, using a phenomenological model. In both weak and strong segregation regimes, large distortions are present in the copolymer film, and these could be observed in experiments.

## II Weakly segregated lamellae

The copolymer order parameter $`\varphi (𝐫)=\varphi _A(𝐫)f`$ is the deviation of the A monomer local volume fraction from its average $`f`$. Above the order-disorder transition (ODT) temperature, the melt is in the disordered, homogeneous state, with $`\varphi (𝐫)=0`$. As the temperature is reduced below the ODT temperature, the system goes through a first order phase-transition to the lamellar phase provided that $`|f\frac{1}{2}|`$ is small enough. Close to the ODT temperature and in the single-mode approximation, the block copolymer (BCP) order parameter is then given by

$$\varphi (𝐫)=\varphi _L\mathrm{cos}(𝐪_\mathrm{𝟎}𝐫)$$
(1)
where $`d_0=2\pi /q_0`$ is the period of lamellar modulations, and $`\varphi _L`$ is their amplitude (to be determined later).

Consider a copolymer melt below the ODT confined by two flat parallel surfaces at $`y=\pm \frac{1}{2}L`$, as in Fig. 1. The surfaces reduce chain entropy, but also chemically interact with the polymers. The difference in the A and B-block surface interactions, $`\sigma _{\mathrm{AS}}`$ and $`\sigma _{\mathrm{BS}}`$, defines the parameter $`\sigma `$,

$$\sigma =\sigma _{\mathrm{AS}}\sigma _{\mathrm{BS}}$$
(2)
In general, $`\sigma (x,y)`$ can be different for the two surfaces, and hence $`\sigma ^+`$ is defined on the $`y=\frac{1}{2}L`$ surface and $`\sigma ^{}`$ is defined for $`y=\frac{1}{2}L`$. The surface interaction, $`F_s`$, can be written as an integral over the bounding surfaces (in units of $`k_BT`$)

$`_s={\displaystyle \left\{\sigma ^{}\varphi ^{}+\sigma ^+\varphi ^+\right\}dxdz}`$ (3)
where $`\varphi ^{}=\varphi (y=\frac{1}{2}L)`$ and $`\varphi ^+=\varphi (y=\frac{1}{2}L)`$. Terms which do not depend on the copolymer order parameter are not important to subsequent calculations, and were dropped out. The interaction in eq 3 is short-range, localized at the surfaces. A positive $`\sigma ^\pm >0`$ induces adsorption of the B monomers ($`\varphi ^\pm <0`$), and $`\sigma ^\pm <0`$ the adsorption of the A monomers ($`\varphi ^\pm >0`$) . We restrict ourselves to homogeneous surfaces, for which $`\sigma ^\pm `$ are taken to be constants over each of the two surfaces.

First let us consider a BCP film without an external electric field. We have already considered this case in a previous publication,$`^{\text{[11]}}`$ and briefly review here the BCP behavior. If the surface affinities $`\sigma ^\pm `$ are sufficiently large, the lamellae will order in a parallel arrangement. These lamellae stretch or compress, increasing the bulk free energy, in order to decrease surface energy. We use below an adaptation of the strong stretching approximation used by Turner $`^{\text{[12]}}`$ and Walton et al. $`^{\text{[13]}}`$ to describe these lamellae. The lamellar period is $`d_0=2\pi /q_0`$ and $`m`$ is the closest integer to $`L/d_0`$. Depending on the values of $`\sigma ^\pm `$, an integer ($`n=m`$) or half integer ($`n=m+\frac{1}{2}`$) number of lamellae exist between the two surfaces. In the former case the ordering is symmetric (the same type of monomers wet both surfaces), while in the latter it is antisymmetric (A monomers wet one surface, whereas B monomers wet the other surface). The parallel lamellae are described by an order parameter $`\varphi _{}`$ given by $`^{\text{[11]}}`$

$$\varphi _{}(y)=\pm \varphi _L\mathrm{cos}[q_{}(y+\frac{1}{2}L)]$$
(4)
The wavenumber is $`q_{}=2\pi n/L`$, and the choice of $`\pm `$ sign in eq 4 is such that the surface interactions, eq 3, are minimized. The amplitude of sinusoidal modulations $`\varphi _L`$ is equal to the amplitude of density modulations in a bulk system, and is given below \[eq 13\].

We consider now the case where an electric field is turned on, in a direction that is perpendicular to the surfaces ($`y`$-axis in Fig. 1). Under conditions of constant voltage difference across the electrodes situated at the two bounding surfaces, the minimum of the free energy is obtained by maximizing the capacitance. Noting that since the A- and B-monomers (blocks) have different dielectric constants, the effect of the electric field is to align the BCP layers parallel to the field, i.e. perpendicular to the surfaces. At a certain field strength, $`E_c`$, this tendency balances the preference for parallel lamellae as induced by the surfaces. Further increase of $`E`$ above the critical value $`E_c`$ gives rise to a perpendicular lamellar ordering.

Close to the ODT, the copolymer ordering is weak, and the energetic cost of compressing or bending the lamellae is small. In this regime the inter-material dividing surface (IMDS) of the lamellar phase, given by the requirement $`\varphi (𝐫)=\frac{1}{2}f`$, can be substantially perturbed from its flat state. This reasoning leads us to the following superposition ansatz. For zero electric field $`E`$, surface interactions orient the lamellae in a parallel orientation. The order parameter is then given by $`\varphi (𝐫)=w(E=0)\varphi _{}(y)`$. The dimensionless amplitude $`w>0`$ is determined by the strength of the surface interactions, and can be larger than unity if $`\sigma ^\pm `$ are sufficiently large. Namely, the surface induced order can be stronger than in the bulk. Upon increase of the electric field, the function $`w(E)`$ of this parallel state diminishes, while the function $`g(E)`$ of perpendicular lamellae (parallel to the electric field) increases. This can be modeled by using the superposition ansatz for the order parameter $`\varphi (𝐫,E)`$ in the presence of the field $`E`$:

$$\varphi (𝐫,E)=w(E)\varphi _{}(𝐫)+g(E)\varphi _{}(𝐫)$$
(5)
The order parameter of the perpendicular lamellae is $`\varphi _{}(x)`$. It depends only on the wavenumber $`q_{}`$, and is given in the single-mode approximation (weak segregation) by

$$\varphi _{}(𝐫)=\varphi _{}(x)=\varphi _L\mathrm{cos}(q_{}x)$$
(6)
The wavenumber $`q_{}`$ is yet to be determined. Without the electric field, $`E=0`$, the amplitude of perpendicular modulations vanishes, $`g=0`$. As the limit $`E\mathrm{}`$ is approached, the effect of the confining surfaces becomes negligible, and the BCP ordering is given by the bulk perpendicular lamellae $`\varphi _{}`$ ($`g=1`$). Hence, the weight amplitudes $`w(E)`$ and $`g(E)`$ satisfy the following limits as function of $`E`$:

$`g(E=0)=0,g(E=\mathrm{})=1`$ (7)
$`w(E=0)=const,w(E=\mathrm{})=0`$ (8)
In order to get explicit expression for the weight functions $`w`$ and $`g`$ we need to consider a specific model. The free energy we use in this section is applicable to the weak segregation regime, and is given by $`=_b+_s`$,$`^{\text{[14, 15, 16]}}`$ with $`_s`$ from eq 3, and $`_b`$ being the bulk contribution. This part of the free energy has a polymer (non-electrostatic) and electrostatic contributions (in units of $`k_BT`$),

$$_b=_p+_{\mathrm{el}}$$
(9)
Hereafter we restrict ourselves to symmetric ($`f=\frac{1}{2}`$) diblock copolymers, where the polymer part of the free energy is $`^{\text{[4, 14, 15, 16, 17, 18]}}`$

$`_p={\displaystyle \left\{\frac{1}{2}\tau \varphi ^2+\frac{1}{2}h\left(q_0^2\varphi +^2\varphi \right)^2+\frac{u}{24}\varphi ^4\right\}\mathrm{d}^3r}`$ (10)
The parameters in eq 10 are given by

$`q_0`$ $`=`$ $`1.95/R_g;h=3\rho c^2R_g^2/2q_0^2`$ (11)
$`\chi _c`$ $``$ $`10.49/N;\tau =2\rho N\left(\chi _c\chi \right)`$ (12)
Denoting $`b`$ as the monomer size, the radius of gyration for Gaussian chains is $`R_g^2\frac{1}{6}Nb^2`$. The polymerization index is $`N`$, the chain density of an incompressible melt is $`\rho =1/Nb^3`$, and $`\chi `$ is the Flory parameter. The amplitude $`\varphi _L`$ in the parallel and perpendicular states is

$$\varphi _L^2=8\tau /u,\tau <0$$
(13)
as is obtained by inserting eq 1 in $`_p`$ of eq 10 and minimizing with respect to $`\varphi _L`$. The dimensionless ratio $`u/\rho `$ and $`c`$ are of order unity, and are taken to be equal exactly to one in the remaining of the paper, $`u/\rho =c=1`$. The correlation length $`\xi 2q_0^2(|\tau |/h)^{1/2}`$ is assumed to be larger than the film thickness, $`\xi L`$.$`^{\text{[19]}}`$

The electrostatic contribution in units of $`k_BT`$ is $`^{\text{[7, 20]}}`$

$`_{\mathrm{el}}=\beta {\displaystyle (\widehat{𝐪}𝐄)^2\varphi _𝐪\varphi _𝐪\mathrm{d}^3q}`$ (14)
$`\beta ={\displaystyle \frac{(\epsilon _A\epsilon _B)^2}{4(2\pi )^4k_BT\epsilon }}`$ (15)
Here $`\varphi _𝐪`$ is the Fourier transform of $`\varphi (𝐫)`$: $`\varphi (𝐫)=\varphi _𝐪\mathrm{exp}(i𝐪𝐫)d𝐪`$, and $`\widehat{𝐪}=𝐪/q`$ is a unit vector in the $`𝐪`$-direction. Copolymer modulations with a non-vanishing component of the wavenumber $`𝐪`$ along the electric field, have a positive contribution to the free energy. In other words, there is a free energy penalty for having dielectric interfaces in a direction perpendicular to the electric field. In eq 15, $`\epsilon _A`$ and $`\epsilon _B`$ are the dielectric constants of the pure A and B-blocks, respectively, $`^{\text{[7, 8]}}`$ and $`\epsilon `$ is the material average dielectric constant for the BCP film,

$$\epsilon =f\epsilon _A+(1f)\epsilon _B$$
(16)
Throughout the remaining of this paper we will focus only on symmetric melts ($`f=\frac{1}{2}`$), having an average dielectric constant $`\epsilon =\frac{1}{2}(\epsilon _A+\epsilon _B)`$. For small concentration variations, valid in the weak segregation regime, $`\epsilon `$ varies linearly with the local copolymer composition $`\varphi `$,

$`\epsilon (\varphi )`$ $`=`$ $`({\displaystyle \frac{1}{2}}+\varphi )\epsilon _A+({\displaystyle \frac{1}{2}}\varphi )\epsilon _B`$ (17)
$`=`$ $`\epsilon +(\epsilon _A\epsilon _B)\varphi `$ (18)
It is possible to perform the spatial integration in eqs. 3, 10 and 14, yielding the free energy per unit volume $`F=F_p+F_{\mathrm{el}}+F_s`$ for a general order parameter $`\varphi (𝐫)`$, eq 5:

$`F`$ $`=`$ $`{\displaystyle \frac{1}{4}}\varphi _L^2\left[\left(\tau +C_{}(E)\right)w^2+\left(\tau +C_{}\right)g^2\right]`$ (19)
$`+`$ $`{\displaystyle \frac{u\varphi _L^4}{64}}(w^4+g^4)+{\displaystyle \frac{u\varphi _L^4}{16}}w^2g^2+{\displaystyle \frac{w\mathrm{\Sigma }}{L}}`$ (20)
The quantities $`C_{}`$ and $`C_{}`$ are positive and given by

$`C_{}(E)`$ $`=`$ $`h(q_0^2q_{}^2)^2+2\beta E^2`$ (21)
$`C_{}`$ $`=`$ $`h(q_0^2q_{}^2)^2`$ (22)
$`\mathrm{\Sigma }=\pm \varphi _L\sigma ^{}\pm \varphi _L\sigma ^+`$ is related to the surface interaction, and is negative. The $`\pm `$ sign is determined from the $`\pm `$ sign of the order parameter in eq 4. The free energies $`F_{}`$ and $`F_{}`$ of the parallel and perpendicular states, respectively, are given as limiting cases

$`F_{}=F(w,g=0)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\varphi _L^2\left[\tau +C_{}(E)\right]w^2+{\displaystyle \frac{u\varphi _L^4}{64}}w^4+{\displaystyle \frac{w\mathrm{\Sigma }}{L}}`$ (23)
$`F_{}=F(w=0,g=1)`$ $`=`$ $`{\displaystyle \frac{1}{4}}\varphi _L^2(\tau +C_{})+{\displaystyle \frac{u\varphi _L^4}{64}}`$ (24)
As discussed in the introduction, we concentrate on the interesting case where in the absence of electric field, the BCP has a parallel ordering given by $`\varphi =w\varphi _{}`$ (with some weight $`w`$). This is equivalent to saying that there exist $`w`$ such that $`F_{}(E=0)<F_{}(E=0)`$. By inserting $`E=0`$ and $`q_{}=q_0`$ in eqs. 23 and 24 we get that in this case

$`{\displaystyle \frac{1}{4}}\varphi _L^2\left[\tau +C_{}(0)\right]w^2+{\displaystyle \frac{u\varphi _L^4}{64}}w^4+{\displaystyle \frac{w\mathrm{\Sigma }}{L}}<{\displaystyle \frac{\tau ^2}{u}}`$ (25)
Therefore, the free energy of the parallel lamellae is expected to be lower than the free energy of the perpendicular lamellae. This assumption is valid for strong enough surface interactions, $`\mathrm{\Sigma }`$. In the next section we proceed to find the weight functions $`w(E)`$ and $`g(E)`$ for any $`E>0`$.

## III Results of the weak segregation model

The free energy (eq 20) is minimized with respect to $`w`$ and $`g`$ to yield two coupled algebraic equations:

$`(\tau +C_{})g\tau g^32\tau w^2g=0`$ (26)
$`\left[\tau +C_{}(E)\right]w\tau w^32\tau g^2w+2\mathrm{\Sigma }/L\varphi _L^2=0`$ (27)
As we will see below, there is a critical value of the electric field, $`E_c`$, which separates between the small and large field behavior. One can estimate $`E_c`$ as the field in which the bulk electrostatic energy $`LE^2`$ (favoring perpendicular lamellae) balances the surface interactions $`\mathrm{\Sigma }`$, and hence $`E_cL^{1/2}`$.

Equations 26 and 27 are analyzed separately for small and large electric fields.

(i) Small electric fields: $`E<E_c`$.

For zero electric field, the lamellae are in their parallel state, namely $`g=0`$. The solution with $`g(E)=0`$ and $`w(E)0`$ \[in eqs. 26 and 27\] corresponds to the minimum of the free energy if $`E`$ is below a certain threshold value $`E_c`$, which is to be determined later. Denoting $`g_<(E)`$ and $`w_<(E)`$ as the weight amplitudes for electric fields $`E<E_c`$, they satisfy the following equations:

$`g_<(E)`$ $`=`$ $`0`$ (28)
$`\left[\tau +C_{}(E)\right]w_<(E)\tau w_<^3(E)+2\mathrm{\Sigma }/L\varphi _L^2`$ $`=`$ $`0`$ (29)
(ii) Large electric fields: $`E>E_c`$.

There is a solution to eqs. 26 and 27 with perpendicular lamellae, namely with a nonzero $`g(E)`$. This solution gives the minimum of the free energy above a critical field, $`E>E_c`$. The weight amplitudes $`g_>(E)`$ and $`w_>(E)`$ are then given by:

$`g_>^2(E)={\displaystyle \frac{\tau +C_{}}{\tau }}2w_>^2(E)0`$ (30)
$`\left[\tau +C_{}(E)2C_{}\right]w_>(E)+3\tau w_>^3(E)+2\mathrm{\Sigma }/L\varphi _L^2=0`$ (31)
The above solution is valid provided that $`w_>`$ is small enough, as given by the inequality

$$w_>^2(E)\frac{1}{2}+\frac{C_{}}{2\tau }\frac{1}{2}$$
(32)
Since $`w(E)`$ is a decreasing function, from eq 22 we see that increasing the electric field $`E`$ from zero, the natural mode of the perpendicular state, $`q_{}=q_0`$, is the first to become critical. Namely, the inequality above is obeyed first when $`q_{}`$ is equal to the bulk mode $`q_0`$, and only later by other $`q`$-modes. For this reason we assume hereafter that $`q_{}=q_0`$, yielding $`C_{}=0`$ and $`g_>^2(E)=12w_>^2(E)`$. Thus, BCP modulations in a direction parallel to the surfaces have the bulk (free) periodicity $`d_0`$.

Equation 31 is a cubic equation for $`w_>`$ and has an analytical solution. It is convenient to express the solution via a parameter $`\theta `$ defined as

$$\mathrm{cos}\left[\theta (E)\right]\frac{\mathrm{\Sigma }/L\varphi _L^2}{3\tau }\left(\frac{\tau C_{}(E)}{9\tau }\right)^{3/2}<0$$
(33)
and take it to be in the range $`\pi <\theta <\frac{3}{2}\pi `$. Recalling that $`\mathrm{\Sigma }`$ and $`\tau `$ are negative, the solution $`w_>`$ to eq 31 is then simply given by

$$w_>(E)=2\sqrt{\frac{\tau C_{}(E)}{9\tau }}\mathrm{cos}\left(\frac{\theta }{3}\right)$$
(34)
Curves of $`g(E)`$ and $`w(E)`$ are shown in Fig. 2 within the assumption $`q_{}=q_0`$. For zero electric field, the film has parallel lamellae, $`g(0)=0`$. With increasing electric field $`E`$, the amplitude $`w_<(E)`$ of parallel lamellae decreases monotonically, as is given by eq 29, while $`g_<(E)`$ remains zero. At the critical field, $`E=E_c`$, there is a first-order phase transition, and perpendicular lamellae appear. The weight amplitudes $`w(E)`$ and $`g(E)`$ are discontinuous, the jump in their values is determined by the degree of segregation $`N\chi `$, inter-surface separation $`L`$ and surface parameters $`\sigma ^\pm `$. Further increase of the electric field causes $`g(E)`$ to increase while $`w(E)`$ decreases. For very large electric fields, $`E\mathrm{}`$, the perpendicular state saturates to its value $`g=1`$, and parallel lamellae completely disappear, $`w=0`$.

The critical field $`E_c`$ is determined by the condition

$$F(w_<,g_<)=F(w_>,g_>)$$
(35)
where the free energy $`F`$ is taken from eq 20. Namely, it is the field where the two values of the free energy cross. The critical field $`E_c`$ as a function of surface separation $`L`$ is shown in Fig. 3 a for weak segregations ($`N\chi =11`$). As the surface separation $`L`$ increases $`E_c`$ decreases with typical oscillations of period $`d_0`$. These oscillations are caused by the frustration occurring when the surface separation is incommensurate with the lamellar period. Figure 3 b shows a log-log plot of the same curve. The dashed line shows a fit to a $`E_cL^{1/2}`$ scaling. This scaling can be obtained by balancing the surface energy $`\mathrm{\Sigma }`$ with the electrostatic contribution $`LE_c^2`$. There is good fit between this $`L^{1/2}`$ line and the peak positions, as expected for unstrained films.$`^{\text{[21]}}`$ However, for frustrated films (where the film thickness is incommensurate with the lamellar period) the deviation from $`E_cL^{1/2}`$ becomes increasingly important as the surface separation $`L`$ is reduced below roughly $`6d_0`$ for the parameters used.

We define $`w_<^c`$ as the value of $`w`$ just before the transition ($`EE_c`$), and $`w_>^c`$ the value after the transition ($`EE_c`$), and similarly for $`g`$. The jump in the weight amplitudes $`\mathrm{\Delta }w`$ and $`\mathrm{\Delta }g`$ is

$`\mathrm{\Delta }w`$ $`=`$ $`|w_>^cw_<^c|`$ (36)
$`\mathrm{\Delta }g`$ $`=`$ $`g_>^cg_<^c=g_>^c`$ (37)
These quantities are plotted as a function of surface separation $`L`$ in Fig. 3 c, for the same parameters as in Fig. 3 a. As $`L`$ increases, the critical field $`E_c`$ decreases, the jump in $`w`$ and $`g`$ gets larger, and the transition from parallel to perpendicular lamellae becomes more abrupt.

In Fig. 4 we show how the film changes its orientation and morphology as the electric field increases. In part a we show a contour plot of the copolymer order parameter $`\varphi =w(E)\varphi _{}+g(E)\varphi _{}`$, for $`E<E_c`$, but only slightly below it. The ordering in the film is parallel to the surfaces, as $`g(E)0`$. In part b the field is slightly increased above its threshold value $`E_c`$, and the undulations created by the appearance of the perpendicular state are prominent. In contrast to the classical Helfrich-Hurault undulations $`^{\text{[22]}}`$ calculated for smectic and cholesteric liquid crystals, here the lateral wavelength is finite, and is equal to the free periodicity $`d_0`$. Note that the correlation length $`\xi `$ is larger than the film thickness, $`\xi >L`$. In addition, the modulations of the adjacent inter-material dividing surfaces (IMDS), given by $`\varphi (𝐫)=0`$, are shown in part c to be out-of-phase with each other and are pronounced only for $`EE_c`$ . Although the perpendicular ordering may be very strong, some parallel ordering is still present (finite $`w>0`$). As the electric field is further increased, there is only little reminiscence of the parallel ordering, and the perpendicular lamellae are nearly perfect. This can be seen in Fig. 4 d, where we choose $`E=4E_c`$.

## IV Strongly segregated lamellae

In this section we consider the same alignment phenomenon under an electric field as in the previous section, but the BCP melt is assumed to be in the strong segregation limit, i.e. $`N\chi N\chi _c10.5`$. In this regime the lamellae are not easily deformed, and the effect of the surface field is important only close to the confining walls, in contrast to the long-range ordering induced by the surfaces in the weak segregation.

### A Unstrained films

We discuss first the alignment phenomenon ignoring the effect of incommensurability between the film thickness $`L`$ and the free lamellar period $`d_0`$. In the following section generalization to strained films will be presented as well. For simplicity, we also assume that the two walls are chemically identical, and preferring the B monomers: $`\sigma _{\mathrm{BS}}<\sigma _{\mathrm{AS}}`$, $`\sigma ^\pm >0`$. For sufficiently strong $`E`$ fields we expect to nucleate a region in the middle of the film with lamellae perpendicular to the walls, as is shown schematically in Fig. 5. Hence two regions of a “T-junction” morphology will exist in the film, in the vicinity of the bounding surfaces. A positive energy penalty (per unit area) $`\gamma __T`$ is associated with each of the two T-junction defects. In principle, other types of defects might exist in the film, but this can only affect the value of $`\gamma __T`$ and not the system behavior as is described below.

We denote the size of the region that is perpendicular to the surfaces by $`l`$. If $`l=0`$ the film has only parallel ordering, while for $`l=L`$ the perpendicular ordering spans the whole film. The free energies per unit area, $`F_{}`$ and $`F_{}`$, for these two extreme cases are:

$`F_{}`$ $`=`$ $`LF_p+2\sigma _{\mathrm{BS}}\lambda _{}LE^2`$ (38)
$`F_{}`$ $`=`$ $`LF_p+\sigma _{\mathrm{AS}}+\sigma _{\mathrm{BS}}\lambda _{}LE^2`$ (39)
$`F_p`$ is the polymer free energy per unit volume of a BCP in the lamellar phase, and in the above we have used two different depolarization factors, $`\lambda _{}`$ and $`\lambda _{}`$, for the two orientations

$`\lambda _{}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{\epsilon _A\epsilon _B}{\epsilon _A+\epsilon _B}}`$ (40)
$`\lambda _{}`$ $`=`$ $`{\displaystyle \frac{1}{16\pi }}(\epsilon _A+\epsilon _B)`$ (41)
These two depolarization factors can be obtained by calculating the electrostatic energy $`(8\pi )^1\epsilon E^2\mathrm{d}^3r`$ of the film and using the boundary condition that the displacement field $`D=\epsilon E`$ is continuous across the two dielectric boundaries. Note that $`\lambda _{}>\lambda _{}`$, meaning that $`\lambda _{}E^2<\lambda _{}E^2`$, and the perpendicular state is favored for $`E\mathrm{}`$.

Let us focus on the mixed state as is illustrated in Fig. 5. For intermediate values of $`l`$, $`0<l<L`$, there are two surface regions of parallel orientation and a central region of perpendicular orientations. The free energy of this mixed state $`F_M`$ (per unit area) is:

$$F_M=LF_p+2\sigma _{\mathrm{BS}}+2\gamma __T\lambda _m(l)LE^2$$
(42)
The system can be regarded as being composed of two capacitors of parallel lamellae and one capacitor of perpendicular lamellae connected in series, yielding the constant $`\lambda _m(l)`$

$$\lambda _m(l)=\left[\frac{Ll}{L\lambda _{}}+\frac{l}{L\lambda _{}}\right]^1$$
(43)
We note that a similar expression for the electrostatic part of the mixed state free energy was previously derived by Pereira and Williams. $`^{\text{[21]}}`$

The value of the surface energies $`\sigma _{\mathrm{AS}}`$, $`\sigma _{\mathrm{BS}}`$ and $`\gamma __T`$ must depend on $`l`$. To see this consider, for example, $`Ld_0`$ and $`l=\frac{1}{2}L`$ . In this case $`\gamma __T`$ is some constant. But as $`lL`$, the value of $`\gamma __T`$ must approach zero, because when $`l=L`$ the T-junction does not exist, and the energy associated with it is zero. We denote by $`a`$ the cutoff length, which is the characteristic width at which $`\gamma __T`$ goes to zero. From the same reason $`\gamma __T`$ must also tend to zero as $`l0`$, with another cutoff length. For simplicity we assume that this cutoff length is also equal to $`a`$. Similarly, the surface interaction energies $`\sigma _{\mathrm{AS}}`$ and $`\sigma _{\mathrm{BS}}`$ tend to zero as $`lL`$, with the same cutoff length. The qualitative forms of these parameters are shown in Fig. 6 a. The smooth monotonic decay to zero at $`l=L`$ and $`l=0`$ is only suggestive. $`\gamma __T(l)`$ has oscillations with period $`\frac{1}{2}d_0`$.

In principle, one should minimize $`F_M`$ with respect to the size of perpendicular domain $`l`$. This generalized (mixed) state would correspond to the parallel state when the minimum is at $`l=0`$, and to the perpendicular state when the minimum is at $`l=L`$. However, we do not know from a molecular description how $`\sigma _{\mathrm{AS}}`$, $`\sigma _{\mathrm{BS}}`$ and $`\gamma __T`$ fall off to zero, and therefore use the approximation shown in Fig. 6 b. $`\sigma _{\mathrm{AS}}(l)=\sigma _{\mathrm{AS}}^0`$ and $`\sigma _{\mathrm{BS}}(l)=\sigma _{\mathrm{BS}}^0`$ are constant for $`l<La`$ and zero for $`l>La`$. $`\gamma __T(l)=\gamma __T^0`$ is constant in the range $`a<l<La`$, and zero otherwise. To make the notation simpler we drop hereafter the superscript zero of $`\sigma _{\mathrm{AS}}`$, $`\sigma _{\mathrm{BS}}`$ and $`\gamma __T`$. Note that the effect of incommensurability between the surface spacing $`L`$ and the lamellar period $`d_0`$, appears as undulations in the plot of $`\gamma __T(l)`$ (Fig. 6 a), and these undulations are neglected here and will be addressed below separately.

With the above assumptions, the size $`l`$ of a perpendicular domain that minimizes the mixed free energy $`F_M`$ is $`l=L2a`$. The free energy of the mixed configuration is thus taken from eq. 42 with $`\lambda _m`$ given by

$$\lambda _m=\left[\frac{2a}{L\lambda _{}}+\frac{L2a}{L\lambda _{}}\right]^1$$
(44)
Using the assumption that the B-polymer is adsorbed at the surfaces in the parallel state ($`\sigma _{\mathrm{BS}}<\sigma _{\mathrm{AS}}`$), we plot in Fig. 7 $`F_{}`$, $`F_{}`$ and $`F_M`$ as a function of electric field strength $`E`$. The dimensionless parameter $`\delta `$ measuring the difference in A and B-block surface interactions is defined as:

$$\delta \frac{\sigma }{\gamma __T}=\frac{\sigma _{\mathrm{AS}}\sigma _{\mathrm{BS}}}{\gamma __T}$$
(45)
Based on the value of $`\delta `$, we now discuss two cases:

(i) Strong surface fields: $`\delta >\delta ^{}`$.

If $`\delta `$ is larger than a threshold value $`\delta ^{}`$, given by

$`\delta ^{}`$ $``$ $`2{\displaystyle \frac{\lambda _{}\lambda _{}}{\lambda _m\lambda _{}}}`$ (46)
$`=`$ $`2+{\displaystyle \frac{4\lambda _{}}{\lambda _{}}}{\displaystyle \frac{a}{L}}+𝒪\left({\displaystyle \frac{a^2}{L^2}}\right)`$ (47)
there are two distinct critical fields $`E_1`$ and $`E_2>E_1`$ given by

$`E_1`$ $`=`$ $`\left[{\displaystyle \frac{2\gamma __T}{L(\lambda _m\lambda _{})}}\right]^{1/2}`$ (48)
$`E_2`$ $`=`$ $`\left[{\displaystyle \frac{\gamma __T(\delta 2)}{L(\lambda _{}\lambda _m)}}\right]^{1/2}`$ (50)
The smallest of these two fields, $`E_1`$, obtained for $`F_M=F_{}`$, (Fig. 7 a) is the field required to create the T-junction defect in the film. This field is independent of $`\delta `$, and for thick films ($`La`$) it scales as $`E_1L^{1/2}`$. For $`E`$ fields in the range $`E_1<E<E_2`$, the film has a region of size $`l=L2a`$ with perpendicular lamellae, while parallel lamellae are localized in a small region of size $`a`$ near the two walls. The second field, $`E_2`$, obtained for $`F_M=F_{}`$, is larger than $`E_1`$ and corresponds to the electric field that is required to destroy the parallel surface layer (of width $`a`$). Note that although these fields are large (typically between $`1`$$`30`$ V/$`\mu `$m), they can be readily achieved in thin-film experiments . For thick films $`E_2`$ is independent of $`L`$ and obtains the asymptotic value

$$E_2\left[\frac{\gamma __T(\delta 2)\lambda _{}}{2a(\lambda _{}\lambda _{})\lambda _{}}\right]^{1/2}$$
(51)
Note also the appearance of $`\lambda _m\lambda _{}`$ and $`\lambda _{}\lambda _m`$ in the denominator of the two critical fields in eqs 48 and 50. If the dielectric contrast between A and B domains is small, namely $`\epsilon _A/\epsilon _B1`$, then $`\lambda _{}\lambda _{}\lambda _m`$, and the critical fields are large. If, on the other hand, $`\epsilon _A/\epsilon _B1`$, then the critical fields required to achieve mixed and perpendicular lamellae are small. For many polymer surface combinations, $`\gamma __T`$ is large and therefore $`\delta `$ is small. In this case the system is classified as having weak surface fields as is discussed below.

(ii) Weak surface fields: $`\delta <\delta ^{}`$.

For small values of $`\delta `$, namely $`\delta <\delta ^{}`$, the free energies $`F_{}`$ and $`F_{}`$ are always smaller than $`F_M`$. Namely, the mixed state is not a possible minimum of the film free energy. There is only one transition at a critical field $`E_3`$ occurring when $`F_{}`$ intersects $`F_{}`$ (Fig. 7 b). This field is given by:

$`E_3`$ $`=`$ $`\left[{\displaystyle \frac{\gamma __T\delta }{L(\lambda _{}\lambda _{})}}\right]^{1/2}`$ (52)
$`=`$ $`\left[{\displaystyle \frac{\sigma }{L(\lambda _{}\lambda _{})}}\right]^{1/2}\delta ^{1/2}L^{1/2}`$ (53)
When $`E=E_3`$, a direct transition occurs from a state where the whole film has parallel lamellae ($`E<E_3`$) to a state where the whole film has perpendicular lamellae without any surface regions ($`E>E_3`$).

The system behavior can be summarized in a phase diagram which depends on the three system parameters: the electric field $`E`$, the surface interaction parameter $`\delta `$ and the film thickness $`L`$. In Fig. 8 a, we show a 2-dimensional cut through the phase diagram, varying $`E`$ and $`\delta `$ while keeping $`L`$ constant. For small electric fields, the lamellae are in the fully parallel configuration. If $`\delta `$ is small, a first-order phase transition to the fully perpendicular state occurs when $`E`$ is increased above $`E_3`$, eq 52. If $`\delta `$ is large enough, namely $`\delta >\delta ^{}`$, there are two transitions when the field is increased. In the regime $`E_2>E>E_1`$, the film is in a mixed state, and layers of thickness $`a`$ with parallel lamellae still exist close to the surfaces. As $`E`$ is further increased above $`E_2`$, the film has a fully perpendicular state. Note that when $`\delta >\delta ^{}`$, $`E_1`$ is the field required to initiate the T-junction defect, and so is independent of $`\sigma `$, while $`E_2`$ is the field that destroys the surface layer, and therefore $`E_2(\sigma 2\gamma __T)^{1/2}`$.

As we have seen, the condition $`\delta >\delta ^{}`$ is required for the existence of the mixed state, and hence for the existence of two critical fields $`E_1`$ and $`E_2`$. Alternatively, this condition can be viewed as a restriction on the film thickness $`L`$ for a given $`\delta `$. Therefore, a mixed state does not exist if

$$L<L^{}\frac{2a(\lambda _{}\lambda _{})[2\lambda _{}+(\delta 2)\lambda _{}]}{\lambda _{}[(\delta +2)\lambda _{}+(\delta 2)\lambda _{}]}$$
(54)
In the limit $`\delta 1`$ and in the case of polystyrene ($`\epsilon 2.5`$) and polymethyl methacrylate ($`\epsilon 6`$) diblock copolymer, one obtains $`L^{}=2a(\lambda _{}\lambda _{})/(\lambda _{}+\lambda _{})1.3a`$. It is smaller for smaller values of $`\delta `$. Hence the mixed state exists for film thickness larger than the range of surface induced ordering $`a`$. We see here again why this mixed state is not expected to occur in the weak segregation limit.

An alternative cut through the phase diagram is shown in Fig. 8 b, where $`\delta `$ is fixed while $`E`$ and $`L`$ are allowed to vary. For small surface separations, $`L<L^{}`$, there is a transition from parallel to perpendicular lamellae at $`E=E_3`$ \[eq 52\]. For larger separations, $`L>L^{}`$, the mixed state appears at $`E=E_1`$ \[eq 48\]. The transition to a fully perpendicular state occurs at $`E=E_2`$ \[eq 50\]. The main difference between this diagram and the one obtained by Pereira and Williams,$`^{\text{[21]}}`$ is that we always find the perpendicular state to be favored for large enough electric fields.

### B Strained films

Strained films occur when the film thickness $`L`$ does not match the lamellar period $`d_0`$, namely $`L/d_0`$ is not an integer or a half-integer number. The effect of this mismatch, which was neglected in the previous section, is considered below. The strong segregation theory of Turner $`^{\text{[12]}}`$ and Walton et al $`^{\text{[13]}}`$ has been successful in describing confined lamellae, and therefore we use it in this section to include the effect of lamellae frustration, and to modify the free energy $`F_{}`$ from section IV A. We use the convention that $`m`$ is the closest integer to $`L/d_0`$, yielding $`F_{}`$ which is the minimum of $`F_{}^s`$ and $`F_{}^{\mathrm{as}}`$, the symmetric and antisymmetric parallel states free energies, respectively:

$`F_{}^s`$ $`=`$ $`\sigma _{\mathrm{AB}}\left[\left({\displaystyle \frac{L}{d_0}}\right)^3{\displaystyle \frac{1}{n^2}}+2n\right]+2\sigma _{\mathrm{BS}}\lambda _{}LE^2`$ (56)
$`\mathrm{for}n=m,\mathrm{symmetric}`$
$`F_{}^{\mathrm{as}}`$ $`=`$ $`\sigma _{\mathrm{AB}}\left[\left({\displaystyle \frac{L}{d_0}}\right)^3{\displaystyle \frac{1}{n^2}}+2n\right]+\sigma _{\mathrm{AS}}+\sigma _{\mathrm{BS}}\lambda _{}LE^2`$ (59)
$`\mathrm{for}n=m\pm {\displaystyle \frac{1}{2}},\mathrm{antisymmetric}`$
Note that in the above expressions for the free energy, eqs 56 and 59, the electrostatic part is exactly like in Sec. IV A. In the symmetric state, both surfaces are wetted by B monomers, while in the antisymmetric state A monomers adsorb to one surface and B monomers adsorb to the second surface.

In Fig. 9 we present the phase diagrams when the mismatch between $`L`$ and $`d_0`$ is taken into account, plotted with same parameters as in Fig. 8. For very thin films, $`L<d_0`$, the perpendicular state is favored over the parallel state, even for $`E=0`$. Apart from this, Fig. 9 a is similar to Fig. 8 a. However, Fig. 9 b is different than Fig. 8 b. The critical fields $`E_1`$ and $`E_3`$ separating the parallel state from the other two states have oscillations with period $`d_0`$. These curves are similar to the critical field curve for weakly segregated lamellae, Fig. 3, computed in Sec. 3, and also to the curve calculated by Ashok et al. using a similar model and different parameters.$`^{\text{[23]}}`$ Note that from the discussion of unstrained films in section IV A, $`E_1`$ and $`E_3`$ have oscillations around lines which decay as $`L^{1/2}`$ for large $`L`$. The field $`E_2`$ separating the mixed state with the perpendicular state depends only weakly on $`L`$.

## V Conclusions

In this paper we study the influence of an applied electric field on the morphology of thin film diblock copolymers. In the absence of an electric field the lamellae are taken to be ordered parallel to the confining surfaces. Strong enough $`E`$ fields in the direction perpendicular to the surfaces will eventually orient the lamellae in a perpendicular direction. However, the response of weakly segregated lamellae is different than the response of strongly segregated lamellae. In the former case, the BCP order parameter is obtained as a function of electric field strength $`E`$ and other system parameters (such as the degree of segregation $`N\chi `$, strength of surface interactions $`\sigma ^\pm `$ and film thickness $`L`$). The field applied in the perpendicular direction diminishes the amplitude of the parallel BCP state. Above the critical field, $`E>E_c`$, a first-order phase transition occurs, from the parallel into the perpendicular state.

The first $`q`$-mode which becomes stable in the perpendicular state $`q_{}`$ has a finite periodicity, equal to the bulk spacing $`d_0=2\pi /q_0`$. Moreover, modulations of adjacent inter-material dividing surfaces (IMDS) are out-of-phase with each other. This is different than the strong segregation instability studied by Onuki and Fukuda,$`^{\text{[20]}}`$ where the slowly varying phase $`\varphi =\varphi _L\mathrm{cos}[q_0x+u(y)]`$ is used, and the free energy is expanded in small $`u`$. At the onset of instability, they found that adjacent IMDS lines are in-phase with each other. The weak segregation IMDS undulations studied here appear because in this regime it is relatively easy to deform the lamellae. For $`E>E_c`$, parallel and perpendicular lamellae coexist in the film. This superposed state is yet to be verified in experiments.

The critical field $`E_c`$ as a function of inter-surface separation $`L`$ decays with characteristic oscillations of period $`d_0`$. These oscillations are the result of lamellar frustration, occurring when the period $`2\pi /q_{}`$ is different than $`d_0`$ (while $`q_{}`$ is taken always to be equal to $`q_0`$). The deviation of the critical field $`E_c`$ from the $`E_cL^{1/2}`$ scaling (of unstrained films) is important only for small surface separations. For large $`L`$ values, the $`L^{1/2}`$ scaling describes well the system behavior.

The weak segregation treatment we present relies on mean-field theory. It is valid close to the critical point, but not too close where critical fluctuations become important.$`^{\text{[24]}}`$ We have ignored the deviations of the inter-material dividing surface from the perfect flat shape that occur near the surfaces. These deviations can be important for small surface separations ($`d_0L`$), or close to the ODT,$`^{\text{[15, 16]}}`$ and should be properly accounted for.

In the second part of the paper we consider BCP film in the strong segregation regime. In contrast to the weak segregation case, here the surface induced ordering has a finite range of length $`a`$. This added length (taken in the weak segregation limit to be much larger than the film thickness $`L`$) results in the existence of two critical fields. Provided that the parameter $`\delta `$ is large enough, $`\delta >\delta ^{}`$ \[eq 46\], for small electric fields, $`E<E_1`$, the film has parallel lamellae. As the field is increased above $`E_1`$ but below $`E_2`$, a region of size $`L2a`$ of perpendicular lamellae nucleates in the middle of the film. In this mixed state, a layer of parallel lamellae of thickness equal to the cutoff length $`a`$ still persists close to the surfaces. When $`E>E_2`$ the system is in the perpendicular state.

The full phase diagram in the $`E`$-$`\delta `$ and the $`E`$-$`L`$ planes is given in Fig. 8. As discussed above, for large $`\delta >\delta ^{}`$ there is a transition from parallel to mixed lamellae at $`E=E_1`$ \[eq 48\] and from mixed to perpendicular lamellae at $`E=E_2`$ \[eq 50\]. Another scenario of a direct transition from a parallel to a perpendicular state at $`E=E_3`$ \[eq 52\] is realized for small $`\delta `$, $`\delta <\delta ^{}`$. These diagrams do not take into account the frustration effects caused by a mismatch between the wall separation and the natural lamellar periodicity in the bulk. In this respect, they are qualitatively applicable to non-polymeric systems such as ferrosmectics in magnetic fields,$`^{\text{[25]}}`$ and to BCP in the hexagonal phase, as investigated experimentally. $`^{\text{[10]}}`$ The difference between the hexagonal and lamellar phase behavior comes from the different value of the depolarization factors $`\lambda _{}`$, $`\lambda _{}`$, and different numerical coefficients in eqs. 38, 39 and 42. These modifications only change the critical fields $`E_1`$, $`E_2`$ and $`E_3`$ by a numerical factor, but do not alter the phase behavior. Indeed, two critical fields between which a regime of mixed state exists have been found experimentally for a BCP is its hexagonal phase by Russell and co-workers.$`^{\text{[10]}}`$ In their experiment, the threshold field (called here $`E_2`$) for transition from the mixed to the fully perpendicular state has no noticeable dependence on surface separation $`L`$, in accord with our calculations as shown in Fig. 8 b and Fig. 9 b.

When the surface separation $`L`$ is not an integer or half integer number of the bulk lamellar period $`d_0`$, the phase diagrams are given by Fig. 9. While in the $`E`$-$`\delta `$ plane the behavior is only slightly changed, in the $`E`$-$`L`$ plane the border of the parallel state has prominent oscillations with period $`d_0`$. Finally, we have not considered the copolymer density undulations as was done by Onuki and Fukuda.$`^{\text{[20]}}`$ These might change the boundary lines between the parallel, perpendicular and mixed phases.

The main prediction of this paper is the dependence of the transition fields on the block copolymer film thickness $`L`$ and the phase diagrams of Fig. 9. It will also be of interest to check experimentally under what conditions the mixed state exists. In addition, the energy penalty $`\gamma __T`$ of the defect created in the film as well as the difference in surface energies of the two blocks, $`\sigma =\sigma _{\mathrm{AS}}\sigma _{\mathrm{BS}}`$, can be deduced in experiment from the measured values of $`E_1`$, $`E_2`$ and $`E_3`$ using eqs. 45, 48 and 50.

## Acknowledgments

We would like to thank T. Taniguchi and T. Yukimura for critical remarks. We benefitted from fruitful discussions with B. Ashok, M. Cloitre, J. DeRouchey, M. Doi, L. Leibler, M. Muthukumar, R. Rosensweig, T. P. Russell, M. Schick, F. Tournilhac and T. Thurn-Albrecht. Partial support from the U.S.-Israel Binational Foundation (B.S.F.) under grant No. 98-00429 and the Israel Science Foundation founded by the Israel Academy of Sciences and Humanities — centers of Excellence Program is gratefully acknowledged.

* Figure 1: Schematic illustration of the system. The two confining surfaces are at $`y=\pm \frac{1}{2}L`$, and have surface interaction parameters $`\sigma ^\pm `$. The electric field points in the perpendicular $`y`$ direction, and is produced by a potential difference between the two surfaces (electrodes).
* Figure 2: Weight functions $`w(E)`$ and $`g(E)`$. The horizontal dash-dot line $`w=1/\sqrt{2}`$ is the value of $`w`$ below which a nonzero $`g(E)`$ is possible. Surface separation is $`L=8d_0`$ and the Flory parameter is $`N\chi =11`$. The film is symmetric, $`\sigma ^+=\sigma ^{}=0.6hq_0^3\varphi _L`$.
* Figure 3: (a) Critical field $`E_c`$ \[in units of $`(hq_0^4/\beta )^{1/2}`$\] required for the appearance of perpendicular lamellae, as a function of surface separation $`L`$. (b) log-log plot of the same curve as in part a. The dashed line is a straight line corresponding to $`E_c0.5(L/d_0)^{1/2}`$ (see text). (c) The jumps $`\mathrm{\Delta }g`$ (dashed line) and $`\mathrm{\Delta }w`$ (solid) in the amplitudes $`g`$ and $`w`$ at the critical field $`E_c`$, as a function of surface separation $`L`$. The film is symmetric, and $`\sigma ^\pm `$ and $`N\chi `$ are as in Fig. 2.
* Figure 4: Contour plots of the BCP order parameter $`\varphi (x,y)=w(E)\varphi _{}(y)+g(E)\varphi _{}(x)`$ for symmetric film. The surfaces are at $`y=\pm \frac{1}{2}L=\pm 2d_0`$, and the field is in the $`y`$ direction. In part a the field is a little smaller than the critical field, $`E=0.98E_c`$, and the film has a perfect parallel ordering. In part b the field is just above the threshold, $`E=1.02E_c`$. The film morphology is a superposition of parallel and perpendicular lamellae. (c) A plot of the IMDS \[given by $`\varphi (x,y)=0`$\] of part b. In part d $`E=4E_c`$, and the lamellae are in the perpendicular state with small distortions. The surface fields are $`\sigma ^+=\sigma ^{}=0.5hq_0^3\varphi _L`$, and the Flory parameter is $`N\chi =11`$ corresponding to correlation length $`\xi 4.8d_0>L`$. The B monomers (colored black) are attracted to the two symmetric surfaces.
* Figure 5: Illustration of the lamellae in the film in the strong segregation regime. A region with perpendicular lamellae exists in the middle of the film, with size $`l`$. Two regions of size $`a`$ each are the parallel lamellar regions.
* Figure 6: (a) Qualitative dependence of $`\gamma __T`$ and $`\sigma _{\mathrm{BS}}`$ on the thickness of perpendicular domain $`\mathrm{}`$, defined in Fig. 5. (b) Simplified curves of part a. The cutoff width is $`a`$, the total film thickness is $`L`$, and $`\gamma __T^0`$ and $`\sigma _{BS}^0`$ are the mean values of $`\gamma __T`$ and $`\sigma _{BS}`$, respectively, for $`a<\mathrm{}<La`$.
* Figure 7: Sketch of the dependence of $`F_{}`$, $`F_{}`$ and $`F_M`$ on the field $`E`$, in solid line, circles and dashed line, respectively. (a) $`\delta >\delta ^{}`$: There are two critical fields, obtained when $`F_M(E_1)=F_{}(E_1)`$ and $`F_M(E_2)=F_{}(E_2)`$. (b) $`\delta <\delta ^{}`$: $`F_M(E)`$ is always larger than $`F_{}(E)`$, and there is only one critical field obtained when $`F_{}(E_3)=F_{}(E_3)`$. Above this field the most stable state is the perpendicular ordering.
* Figure 8: (a) Phase diagram in the $`E`$-$`\delta `$ plane. If $`\delta =(\sigma _{\mathrm{AS}}\sigma _{\mathrm{BS}})/\gamma __T<\delta ^{}`$, there is a transition between parallel and perpendicular lamellae at $`E=E_3`$. For $`\delta >\delta ^{}`$, there is a transition from the parallel to the mixed state at $`E=E_1`$, followed by a second transition from the mixed to the perpendicular state when $`E=E_2>E_1`$. Surface separation is chosen as $`L=10d_0`$ and $`a=d_0`$. (b) Similar diagram, but in the $`E`$-$`L`$ plane, with $`\delta =5`$. In both parts, the effect of mismatch between $`L`$ and $`d_0`$ is neglected, and the electric fields are scaled by $`(\gamma __T/d_0)^{1/2}`$.
* Figure 9: Phase diagrams as in Fig. 8, but for strained films. The critical fields $`E_1`$ and $`E_3`$ have oscillations as a function of $`L`$, in contrast to Fig. 8 b, but follow on average a similar $`L^{1/2}`$ scaling behavior.


format: text

text:
# V2487 Ophiuchi: A strong candidate for a recurrent nova and a progenitor of Type Ia supernova

## 1. Light Curve Analysis and a Very Massive White Dwarf

The outburst of V2487 Ophiuchi (Nova Oph 1998) was discovered at $`m_V9.5`$ by Takamizawa on 1998 June 15.561 UT. The early optical decline rate of 0.37 mag day-1 makes this nova one of the fastest ever seen (see Fig. 1, from the VSNET archives). The rapid decline stopped about 10 days after maximum and then the decline rate lessened to $`0.05`$ mag day-1, that is, the brightness stayed at $`m_V14`$ from 10 to 30 days after maximum. We call this period mid-plateau phase. Then the optical brightness again declined so fast from $`m_V14`$ to $`m_V16`$ in about 10 days. The brightness stayed again at $`m_V16`$ at least until 70 days after maximum (HJD 2,451,050), which is not shown in Fig.1 but from observational points reported in the AAVSO archives.

We have modeled the system consisting of a very massive WD and a lobe-filling main-sequence (MS) star. Irradiation effects of the accretion disk (ACDK) and the MS companion by the WD are included into the light curve calculation. The numerical method has been described in Hachisu & Kato (2001). We are able to reproduce the light curve by adopting model parameters similar to those of U Sco (Hachisu et al. 2000). The model parameters are shown in Fig. 1.

The very rapid decline and the ensuing mid-plateau phases are common features among the U Sco subclass of recurrent novae. If V2487 Oph belongs to the U Sco subclass, its orbital period should be between $`0.3`$ and $`3`$ days. The numerical results strongly indicate that this nova is a recurrent novae. Adding V2487 Oph to the member of the recurrent novae would provide us valuable information on recurrent novae and SN Ia progenitors.

## References

Hachisu, I., & Kato, M. 2001, ApJ, 558, 323

Hachisu, I., Kato, M., Kato, T, & Matsumoto, K. 2000, ApJ, 528, L97


text:
# Estimation of Initial Conditions and Secure Communication

## Abstract

We estimate the initial conditions of a multivariable dynamical system from a scalar signal, using a modified Newton-Raphson method incorporating the time evolution. We can estimate initial conditions of periodic and chaotic systems and the required length of scalar signal is very small. We also find that the information flow from one variable to the other has logarithmic dependence on time. An important application of the method is in secure communications. The communication procedure has several advantages as compared to others using dynamical systems.

A trajectory of a given dynamical system in its state space depends on the set of initial conditions with which it starts. In particular, the state of a chaotic system at a later time is exponentially sensitive to changes in its initial state . This feature of a chaotic system leads to a complex behaviour in state space that appears random yet is deterministic, and the time evolution is uniquely fixed by the initial state of the system. Though there are several invariant measures of a chaotic system which are not sensitive to the initial conditions, the exact trajectory crucially depends on the initial state and hence is difficult to reproduce due to sensitivity to initial conditions.

In light of these facts, it is interesting and important to ask whether the complete set of initial conditions of a given multivariable dynamical system can be estimated from a given scalar time series of a single state space variable. In this letter we present a novel and simple method to estimate the initial conditions from a given scalar time series. The method is based on a modified multidimensional Newton-Raphson method  that includes the time evolution of the system. The length of the time series required for estimating the initial conditions is very small. Also, the method works even when the conditional Lyapunov exponents are positive.

An important application of our method is in the area of secure communications . A nice feature of this application is that the signal that is transmitted is not the one modulated by the information signal, making it difficult to crack the method. Our method is also useful in the problem of synchronization of chaotic signals . The knowledge of initial conditions means the response system can be be synchronized with the drive system almost instantaneously thereby removing the problem of transients. Also synchronization can be achieved in most of the cases where other methods fail .

Let us consider an autonomous dynamical system given by,

$$\dot{𝐱}=𝐅(𝐱),$$
(1)
where $`𝐱=(x_1,x_2,\mathrm{},x_d)`$ is a $`d`$-dimensional state vector whose evolution is governed by the function $`𝐅=(F_1,F_2,\mathrm{},F_d)`$. Given an initial state vector $`𝐱(0)`$ at time $`t=0`$, the time evolution $`𝐱(t)`$ is uniquely determined by Eq. (1). Now let us assume that only one component of the state vector is known to us and we take it to be $`x_1(t)`$ without loss of generality. The problem that we address is to obtain the initial state vector $`𝐱(0)`$ from the knowledge of the scalar signal $`x_1(t)`$.

Let $`𝐲(0)`$ denote a random initial state vector and $`𝐲(t)`$ its time evolution obtained from Eq. (1). Let $`𝐰(t)`$ denote the difference $`𝐰(t)=𝐲(t)𝐱(t)`$. We look for the solution of the equation

$$𝐰(t)=0.$$
(2)
Noting that the initial state vectors $`𝐲(0)`$ and $`𝐱(0)`$ uniquely determine the difference $`𝐰(t)`$, one of the solutions of Eq. (2) is $`𝐲(\mathrm{𝟎})𝐱(0)=0`$ and this is the solution that we are searching for.

We now introduce the notation $`𝐰^n=𝐰^n(𝐲^0,𝐱^0)=𝐰(n\mathrm{\Delta }t)`$, where $`\mathrm{\Delta }t`$ is a small time interval. Similarly, $`𝐲^n=𝐲(n\mathrm{\Delta }t)`$ and $`𝐱^n=𝐱(n\mathrm{\Delta }t)`$. With this notation condition (2) can be written as $`𝐰^n=0`$.

Our approach to the solution of Eq. (2) is a modified Newton-Raphson method  which includes the time evolution of the system.

Let us first consider $`𝐰^1`$. We have

$`0`$ $`=`$ $`𝐰^1(𝐱^0,𝐱^0)=𝐰^1(𝐲^0+\delta 𝐲^0,𝐱^0),`$ (3)
$`=`$ $`𝐰^1(𝐲^0,𝐱^0)+(\delta 𝐲^0\mathbf{}\mathbf{}_{𝐲^\mathrm{𝟎}})𝐰^1(𝐲^0,𝐱^0)+𝒪((\delta 𝐲^0)^2),`$ (4)
where $`\delta 𝐲^0=𝐱^0𝐲^0=𝐰^0`$ and the last step is a Taylor series expansion in $`\delta 𝐲^0`$. For small $`\mathrm{\Delta }t`$, we can write

$$𝐰^1(𝐲^0,𝐱^0)=𝐰^0+\mathrm{\Delta }t[𝐅(𝐲^0)𝐅(𝐱^0)]+𝒪((\mathrm{\Delta }t)^2).$$
(5)
Substituting Eq. (5) in Eq. (4) and neglecting higher order terms, we get,

$$𝐰^1(𝐲^0,𝐱^0)=𝐰^0\mathrm{\Delta }t(𝐰^0\mathbf{}\mathbf{}_{𝐲^\mathrm{𝟎}})𝐅(𝐲^0).$$
(6)
It is convenient to write the above equation in a matrix form as

$`W^1`$ $`=`$ $`(I+\mathrm{\Delta }tJ^0)W^0=A^0W^0,`$ (7)
where $`W^n`$ is the column matrix corresponding to the vector $`𝐰^n`$, $`I`$ is the identity matrix, $`A^n=I+\mathrm{\Delta }tJ^n`$, and the elements of the Jacobian matrix $`J^n`$ are $`J_{ij}^n=F_i(𝐲^n)/y_j^n`$.

Next we consider $`𝐰^2`$ or $`W^2`$. Proceeding as above, we get (see Eq. (7)),

$`W^2`$ $`=`$ $`(I+\mathrm{\Delta }tJ^1)W^1=(I+\mathrm{\Delta }tJ^1)(I+\mathrm{\Delta }tJ^0)W^0`$ (8)
$`=`$ $`A^1A^0W^0.`$ (9)
Similarly, the equation for $`W^n`$ is

$`W^n`$ $`=`$ $`A^{n1}A^{n2}\mathrm{}A^0W^0.`$ (10)
We now concentrate on the first component of the signal whose time series is assumed to be known. For a $`d`$-dimensional system we need $`d1`$ equations to determine the initial state vector $`𝐱^0`$. Eqs. (7), (9) and (10) give us the required relations.

$`W_1^1`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{d}{}}}A_{1i}^0W_i^0,`$ (11)
$`W_1^2`$ $`=`$ $`{\displaystyle \underset{i,j=1}{\overset{d}{}}}A_{1i}^1A_{ij}^0W_j^0,`$ (12)
$`\mathrm{}`$ (13)
$`W_1^{d1}`$ $`=`$ $`{\displaystyle \underset{i,\mathrm{},l,m=1}{\overset{d}{}}}A_{1i}^{d2}\mathrm{}A_{lm}^0W_m^0,`$ (14)
These are $`d1`$ simultaneous equations for $`W^0`$.

The numerical procedure is as follows. We set the initial state of system (1) to a random initial guess vector $`\left(𝐲^0\right)_{old}`$ with $`\left(y_1^0\right)_{old}=x_1^0`$ and evolve it using Eq. (1). Using this vector $`𝐲(t)`$ we write down $`d1`$ simultaneous equations (Eqs. (14)) which can be solved for $`d1`$ unknown components of $`𝐰^0=\delta 𝐲^0`$. Also, $`\delta y_1^0=0`$. Thus the initial guess vector can be improved by

$$\left(𝐲^0\right)_{new}=\left(𝐲^0\right)_{old}+\delta 𝐲^0.$$
(15)
This sets up an iterative scheme giving us better and better estimates of the initial vector which converge to $`𝐱^0`$.

We note that as in Newton-Raphson method, the choice of the initial guess vector can be very important . In some cases, the iterative procedure of Eq. (15) may not converge or converge to a wrong root. In such cases, a different choice of initial guess vector can be useful.

We further note the similarity of our method with the so called method of variational equations in analytical dynamics . The method of variational equations can be applied to a known Hamiltonian system to determine an unknown neighboring trajectory to an already known one. There, the method requires a complete particular solution of a known set of Hamiltonian equations of motion. In contrast, we have used our method for dissipative chaotic systems. In such systems an analytical solution of the equations of motion cannot be known. Further our method requires only one component of a complete trajectory to be sampled. This is important for the application of our method to secure communications as we will demonstrate afterwards.

We now illustrate our method of estimating the initial state. As an example we discuss the Rössler system given by , $`\dot{𝐱}=(x_2x_3,x_1+ax_2,b+x_3(x_1c))`$. First, we consider a case when the time series for $`x_1`$ is given and we want to estimate $`(x_2^0,x_3^0)`$. We chose the parameters $`(a,b,c)`$ such that the system is in the chaotic regime and the initial state $`𝐱^0`$ is on the chaotic attractor. We start with an arbitrary initial state $`𝐲^0=(y_1^0,y_2^0,y_3^0)`$ with $`y_1^0=x_1^0`$. From Eqs. (14) we get a pair of simultaneous equations for $`(\delta y_2^0,\delta y_3^0)`$, which can be solved to obtain these values. With $`\delta y_1^0=0`$ we use these in an iterative manner (Eq. (15)) to obtain the correct initial conditions.

Let $`e_i=|y_i^0x_i^0|`$ denote the absolute error in the estimation of $`x_i^0`$. In Fig. 1(a) we plot a graph of errors $`e_2`$ and $`e_3`$ plotted against the number of iterations, $`n`$, of our method (Eq. (15)). From Fig. 1(a) we see that the errors go to zero and the successive estimates converge to the correct values of $`(x_2^0,x_3^0)`$. Using only two data points in the given time series $`x_1(t)`$, we can thus readily estimate the full initial state $`𝐱^0`$. We also note that the rate of convergence is very good. In about 8 to 10 iterates we obtain the initial values $`(x_2^0,x_3^0)`$ to within computer accuracy. If we write the deviations of the successive iterates from the correct values in the form

$$(e_i)_n=\left|\left(y_i^0\right)_nx_i^0\right|e^{\alpha n},$$
(16)
where $`n`$ is the number of iterations, then the value of the parameter $`\alpha `$ is found to be $`2.12`$ for $`e_2`$ and $`2.10`$ for $`e_3`$. This is consistent with the fact that Newton-Raphson method has a quadratic convergence .

We note that the largest Lyapunov exponent for the subsystem $`(y_2^0,y_3^0)`$ (conditional or subsystem Lyapunov exponent) is positive . The success of our method does not depend on whether this Lyapunov exponent is positive or negative. This is important for synchronization of chaotic signals.

We next present cases where time series for the variables $`x_2`$ and $`x_3`$ of the Rössler system are given. The procedure is similar to the case of time series for $`x_1`$ as discussed above. Fig. 1 (b) shows the errors $`e_1`$ and $`e_3`$, when time series for $`x_2`$ is given, plotted against the number of iterations. The parameter $`\alpha `$ (Eq. (16)) is $`2.08`$ for $`e_1`$ and $`2.10`$ for $`e_3`$. This again indicates a quadratic convergence. Similarly, Fig. 1 (c) shows the quantities $`e_1`$ and $`e_2`$ when time series for $`x_3`$ is given, as a function of the number of iterations. The parameter $`\alpha `$ (Eq. (16)) is $`2.02`$ for $`e_1`$ and 2.10 for $`e_2`$. We note that the largest subsystem Lyapunov exponent is negative when time series for $`x_2`$ is given and is positive when time series for $`x_3`$ is given .

We have successfully applied our method to estimate the initial state vector using a given scalar time series for many other dynamical systems as well. These include Lorenz system  in its periodic, chaotic or intermittent regimes, Chua’s circuit , the disk dynamo system modeling a periodic reversal of earth’s magnetic field , a 3-d plasma system formed by a three wave resonant coupling equations  and a four dimensional phase converter circuit .

We now turn to the important question of the rate of flow of information about other variables into the given variable, i.e. the variable whose time series is known, say $`x_1`$. To investigate this question we start with the initial value $`x_1(0)`$ at $`t=0`$. This alone does not tell us anything about the other variables. As time evolves information about other variables flows into $`x_1`$ and the method presented above allows us to extract this information. As a specific example consider the Rössler system with time series for $`x_1`$. Fig. 2 shows the plot of the asymptotic error $`e_2`$ of the estimation of $`x_2(0)`$ for large $`n`$, denoted by $`e_{\mathrm{}}`$, as a function of the total time interval $`\tau =(d1)\mathrm{\Delta }t`$ for which the time series of $`x_1`$ is used in the calculations. We observe that the accuracy of estimation improves with $`\tau `$ and shows a power law behaviour of the type

$$e_{\mathrm{}}=\tau ^\nu $$
(17)
Fig. 2 shows a very good straight line behaviour on a log–log plot over seven orders of magnitude and the slope gives the value of the exponent $`\nu =2.00\pm 0.01`$. The power law dependence of the accuracy on time interval as in Eq. (17) shows that the information entropy in terms of bits or partitions grows logarithmically with time. Note that this flow of information is different from the usual linear growth of information entropy with time in chaotic systems which corresponds to an exponential dependence of $`e_{\mathrm{}}`$ on $`\tau `$. The excellent straight line fit over seven orders of magnitude that we observe in Fig. 2, makes us believe that the power law dependence of Eq. (17) and the corresponding logarithmic growth of information content with time is a genuine property of dynamical systems. Further, the exponent $`\nu `$ appears to be universal and has the same value within numerical errors for the various examples that we have considered.

Now we will discuss an important application of our method for secure communications. Let us suppose that Eq. (1) describes a $`d`$-dimensional chaotic system at the transmitter end and a replica of the same at the receiver end. Let $`s(i),i=1,2,\mathrm{},k`$ be the information signal to be encoded for communication and without loss of generality let $`x_1(t)`$ be the scalar output signal to be transmitted. The signal $`s(t)`$ is encoded as follows. While choosing the initial conditions of the vector $`𝐱(0)`$, the component $`x_1(0)`$ is chosen arbitrarily while the other $`d1`$ components are chosen as $`x_i(0)=s(i1)`$, $`i=2,\mathrm{},d`$. Using the initial vector, $`𝐱(0)`$, we evolve Eq. (1) for $`d1`$ time steps of $`dt`$ to obtain the first $`d`$ values of the scalar variable $`x_1`$. This is the first stage of evolution. The vector $`𝐱(t_1)`$, where $`t_n=n(d1)dt`$, is modified using $`s(i)`$ as explained below (Eq. (19)) and the resultant is evolved for $`d1`$ time steps completing the second stage. The process is repeated as often as required. After completion of each stage the vector $`𝐱(𝐭_𝐧)`$ is modified as follows.

$`x_1(t_n)`$ $`=`$ $`x_1(t_n)`$ (18)
$`x_i(t_n)`$ $`=`$ $`x_i(t_n)+s(n(d1)+i1),i=2,\mathrm{},d.`$ (19)
Thus, assuming that $`k`$, the number of elements of information signal $`s(i)`$, is divisible by $`d1`$ we get $`k+1`$ values of the encoded signal $`x_1`$.

The recovery of the information signal at the receiver end using our procedure of initial condition estimation is straightforward. At each stage $`d1`$ values of $`x_1`$ unable us to estimate the initial conditions thereby recovering the signal from Eq. (19).

The communication procedure described above has several advantages. Note that (a) The transmitted signal $`x_1`$ is not directly modulated by the information signal $`s`$. (b) The information signal is added to more than one variables. Due to these reasons, local approximation for the flow of the dynamical system in the embedded phase space  cannot be used to attack the coding procedure. Secondly, the transmitted signal is of the same size as that of the information signal, and there is no initial transient as required in other synchronization based methods . Thirdly, as in other methods using chaotic encoding the same signal $`s`$ will give different outputs for reruns of the coding procedure thereby making any attack very difficult.

Another important application of our initial condition estimation method is in the problem of synchronization of two identical chaotic systems coupled unidirectionally by a scalar signal. The knowledge of initial conditions enables us to obtain immediate synchronization thereby eliminating transients . Also, the method could be used repeatedly to maintain synchronization. This will be particularly useful for systems where the largest conditional Lyapunov exponent is positive and other methods of synchronization fail .

To summarize, we have introduced a novel yet simple method to estimate initial conditions of a multivariable dynamical system from a given scalar signal. Our method is based on a multidimensional Newton-Raphson method where we include the time evolution of the system. The method gives a reasonably fast convergence to the correct initial state. The required length of the time series is very small. The method works even when the largest conditional Lyapunov exponent is positive. We find that the information flow between different variables grows logarithmically with time. This is different from the usual linear growth of information in chaotic systems. An important application of the method is in secure communications. This method of communication has several advantages as compared to other methods of communication based on nonlinear dynamical systems.


text:
# AN EMPIRICAL ANALYSIS OF MEDIUM-TERM INTEREST RATES

## 1 Introduction

In financial theory and practice, interest rates are a very important subject which can be approached from several different perspectives.

The classical theoretical approach models the term structure of interest rates using stochastic processes and the no arbitrage argument. Various models have been proposed and can be found in Pagan (1997), Bouchaud (1997), Rebonato (1998), Wilmott (1999), Hull (2000), Mantegna (2000), Dacorogna (2001), Vasicek (1977), Cox (1985), Heath (1992), Hull (1993) . Although they provide analytical formulas for the pricing of interest rate derivatives, the implied deformations of the term structure have a Brownian motion component and are often rejected by empirical data (see Chan (1992) ). The inadequacies of the Gaussian model for the description of financial time series has been reported since a long time ago by Mandelbrot (1963), see also Mandelbrot (1997) , but thanks to the availability of large sets of financial data, the interest on this point has risen recently in Mantegna (2000), Dacorogna (2001), Sornette (2000), Muller (1995), Mantegna (1994), Ghashghaie (1996), Mantegna (1997) . In particular, the fat-tail property of the empirical distribution of price changes has been widely documented and is a crucial feature for monitoring the extreme risks (Embrechts (1997), Embrechts (1999), McNeil (1998), McNeil (1999) ). Many of the recent studies concern high frequency data of stock indexes or exchange rates (Muller (1995), Ghashghaie (1996), Mantegna (1997), Muller (1996) ).

Also for interest rates, up to now, no universally accepted theory has been obtained for the description of experimental data. Dacorogna has investigated the intraday behavior of the interest rate markets (Muller (1996), Piccinato (1997), Ballocchi (1997), Ballocchi (1999) ). Bouchaud et al. (Bouchaud (1997), Matacz (1999), Matacz (2000) ) have considered the interest rate curve as a vibrating string subject to random shocks along its profile. Along these lines, Sornette (1998) , has derived a general condition that the partial differential equations governing the motion of such string must obey in order to account for the condition of absence of arbitrage. Jean Nuyts (1999)  has developed a phenomenological approach based on Padé Approximants. In the work of Jamshidian (1997) , a self-contained theory is presented for pricing and hedging LIBOR and swap derivatives by arbitrage. Appropriate payoff homogeneity and measurability conditions are identified which guarantee that a given payoff can be attained by a self-financing trading strategy. An important recent development in the pricing of interest rate derivatives is the emergence of models that incorporate lognormal volatilities for forward LIBOR or forward swap rates while keeping interest rates stable. Zhao et al. (1999)  have introduced methods for discretizing these models giving particular attention to precluding arbitrage among bonds and to keeping interest rates positive even after discretization.

In this framework, we have empirically studied the probability density distribution of LIBOR, in order to characterize the stochastic behavior of the daily fluctuations.

In Section 2, we present the data set. The data are then analyzed and the results are given in Section 3. Section 4 contains a discussion and some final remarks.

## 2 The data set

LIBOR stands for the London Interbank Offered Rate and is the rate of interest at which banks are willing to offer deposits to other prime banks, in marketable size, in the London interbank market.

BBA (British Bankers’ Association) LIBOR  is the most widely used benchmark or reference rate. It is defined as the interest rate at which deposits are perceived to be generally available in the London interbank at $`11AM`$ (Greenwich time).

Moreover, it is used as the basis for settlement of interest rate contracts on many of the world’s major future and option exchanges as well as most Over the Counter (OTC) and lending transactions.

The BBA, advised by senior market practitioners, maintains a reference panel of $`16`$ banks. The purpose is to produce a reference panel of banks reflecting the balance of the market by country and by type of institution. Individual banks are selected within this principle on the basis of reputation, scale of market activity and perceived expertise in the currency concerned. The BBA surveys the panel’s market activity and publishes their market quotes on-screen. The highest four and lowest four market quotes are disregarded and the remaining eight quotes are averaged: the resulting spot fixing is the BBA LIBOR rate. The quotes from all $`16`$ panel banks are published on-screen to ensure transparency. BBA LIBOR is compiled each working day by an electronic vendor (currently Dow Jones Telerate) and broadcast through ten international distribution networks including Reuters, Dow-Jones Telerate, Knight-Ridder, and Bloomberg.

BBA LIBOR fixings are provided in seven international currencies: Pound Sterling, US Dollar, Japanese Yen, Swiss Franc, Canadian Dollar, Australian Dollar, EURO. LIBOR rates are fixed for each currency at monthly maturities from one month to $`12`$ months. Rates shall be contributed in decimal to at least two decimal places but no more than five.

In the following, we have analyzed a data set of LIBOR interest rates $`r(T,t)`$, where $`T`$ is the maturity date and $`t`$ the current date, for EURO and Pound Sterling. These data are shown in Fig.1 and Fig.2 where $`t`$ goes from January 2, $`1997`$ to September $`17`$, $`1999`$, and $`T`$ assumes the following values: $`1`$, $`3`$, $`6`$, $`9`$ and $`12`$ months for the Pound Sterling and $`1`$, $`3`$, $`6`$, $`12`$ months for the EURO.

In Figs.3, 4 the 1-month LIBOR is compared to the interest rates fixed by Central Banks at that time, namely the REPO (repurchase agreement) data. BBA LIBOR follows the trend determined by the decisions of Central Banks. In order to roughly eliminate these trends, in Figs.5 and 6 the interest rates differences $`\mathrm{\Delta }r(T,t)`$=$`r(T,t+\mathrm{\Delta }t)r(T,t)`$, with $`\mathrm{\Delta }t`$ being $`1`$ day and $`T`$=$`1`$ month, are plotted as a function of the current date, for the EURO and the Pound Sterling, respectively. A similar behavior is also found for the other maturities. Some large oscillations of $`\mathrm{\Delta }r`$ are induced by Central Banks. In any case, $`\mathrm{\Delta }r`$ heavily fluctuates around zero.

## 3 Results

In this section, we focus the attention on the probability distribution behavior of the interest rates increments $`\mathrm{\Delta }r(T,t)`$ defined in the previous section. To this purpose, we estimate $`\mathrm{\Psi }(\mathrm{\Delta }r)`$, the complementary cumulative distribution function of the daily interest rates increments, defined as:

$$\mathrm{\Psi }(\mathrm{\Delta }r)=1_{\mathrm{}}^{\mathrm{\Delta }r}p(\eta )𝑑\eta $$
(1)
where $`p`$ is the probability of $`\mathrm{\Delta }r(T,t)`$.

Because LIBOR data are supplied with only few decimal digits, it is interesting to examine the effects of different data cut-offs in the behavior of $`\mathrm{\Psi }(\mathrm{\Delta }r)`$. In Fig.7, we plot the complementary cumulative distribution function for a simulated Gaussian stochastic process using data characterized by three different decimal digit precisions. It turns out that the numerical rounding does not influence the results.

Figs.8 and 9 show the tail distribution behaviors in the case of EURO and Sterling Pound, respectively. In particular, in Fig.8, we report the empirical results obtained estimating the probability density function of both positive and negative LIBOR increments with $`\mathrm{\Delta }t`$=$`1`$ day and $`T`$=$`1`$ month. These empirical curves are slightly asymmetric and the negative variations are more probable than the positive one. In the same figure, these two curves are compared with the equivalent (i.e. with the same average and standard deviation) Gaussian complementary cumulative distribution. The non-Gaussian behavior is also evident from Fig.9. In both Figs.8 and 9, the empirical LIBOR data exhibit a fat tail or leptokurtic character, which is present for the other maturities as well. These observations indicate that the random behavior of $`\mathrm{\Delta }r(T,t)`$ is non-Gaussian and that using a Gaussian probability density function leads to underestimating the probability of large fluctuations.

For a better understanding of what kind of stochastic process we are dealing with, we present some results on the power spectral density behavior (Kay (1981) ).

The power spectra, $`S(f)`$, for both $`r(1,t)`$ and $`\mathrm{\Delta }r(1,t)`$ are reported in Figs.10 and 11 for the EURO and in Figs.12 and 13 for the Sterling Pound. For $`r`$ the spectral density shows a power law behavior. A linear fit gives a slope value $`\alpha =1.80\pm 0.02`$ and $`\alpha =1.79\pm 0.01`$ for EURO and Sterling Pound, respectively. A similar result holds for the other maturities. Therefore, we argue that the power spectrum analysis for $`r(T,t)`$ indicates a stochastic process with spectral components decreasing as $`S(f)f^\alpha `$ (Feller (1971) ). As for the increments, in Figs.11 and 13, the power spectrum is flat, typical of a white noise process without seasonality peaks. These results are also corroborated by a similar analysis performed on Eurodollars interest rates for a longer time period (Di Matteo (2001) ).

## 4 Discussion

In the previous section, we have shown that the daily increment of the interest rate series is non-Gaussian and follows a leptokurtic distribution. The problem arises of how derivatives written on interest rates can be evaluated, given that the usual Gaussian white-noise assumption of many models is not satisfied. Indeed, a first partial answer, is that the central limit theorem ensures that, after a sufficiently long time, the increment distribution will tend to a Gaussian distribution. However, if the time horizon of derivative evaluation is not appropriate, deviations from the Gaussian behavior may lead to a dramatic underestimate of large increments with a consequent improper risk coverage.

Although well studied in mainstream finance, this problem has received much attention in recent times, within the community of physicists working on financial problems (Bouchaud (1997), Bouchaud (1994) ). In particular, Bouchaud and Sornette (Bouchaud (1994) ) suggested the direct use of the historical probability measure, rather than the equivalent martingale measure for evaluating options. In this way, one gets an option price depending on the expected rate of returns, a consequence which is not fully desirable due to the subjective character of that rate. Assessing trends is a difficult task, as they depend on decision taken by Central Banks (as shown in Figs.3 and 4) and are based on macroeconomic effects. Thus, martingale methods could prove more reliable. As early as 1977, some years after the seminal paper of Black and Scholes, Parkinson generalized their approach to option pricing and explicitly took into account leptokurtic distributions (Parkinson (1977) ). More recently, Boyarchenko and Levendorskii have studied the problem of option pricing in the presence of a specific distribution which seems to fit well the empirical data in many instances: the truncated Lévy distribution (Boyarchenko (2000), Koponen (1995) ). Along these lines, we believe, it will be possible to develop a consistent option pricing model taking into account the leptokurtic character of the empirical short to mid-term interest rate distributions.

## Acknowledgements

We thank the Risk Management division of Banca Intesa (Milano) for providing us the data. T. Di Matteo wishes to thank Prof. Sandro Pace and Tomaso Aste for fruitful discussions and support. This work was partially supported by INFM (Istituto Nazionale per la Fisica della Materia) FRA project - $`AMF_2`$: Application of Physical Methods to Financial Market Analysis.


text:
# Multi Module Model for Ultra-Relativistic Heavy Ion Collisions

## Abstract

The Multi Module Model for Ultra-Relativistic Heavy Ion Collisions at RHIC and LHC energies is presented. It uses the Effective String Rope Model for the calculation of the initial stages of the reaction; the output of this model is used as the initial state for the subsequent one-fluid calculations. It is shown that such an initial state leads to the creation of the third flow component. The hydrodynamical evolution of the energy density distribution is also presented.

A realistic and detailed description of an energetic heavy ion reaction requires a Multi Module Model, wherein the different stages of the reaction are each described with a suitable theoretical approach. It is important that these Modules are coupled to each other correctly: on the interface, which is a three dimensional hyper-surface in space-time, all conservation laws should be satisfied and entropy should not decrease. These matching conditions were worked out and studied for the correct matching at FO in detail in Refs. .

Fluid dynamical models are widely used to describe ultra-relativistic heavy ion collisions. Their advantage is that one can vary flexibly the Equation of State (EoS) of the matter and test its consequences on the reaction dynamics and the outcome. For example, the only models which may handle the supercooled QGP are hydrodynamical models with a corresponding EoS. In energetic collisions of large heavy ions, especially if a Quark-Gluon Plasma (QGP) is formed in the collision, one-fluid dynamics is a valid and good description for the intermediate stages of the reaction. Here, interactions are strong and frequent, so that other models, (e.g. transport models, string models, etc., that assume binary collisions, with free propagation of constituents between collisions) have limited validity.

On the other hand, the initial and final, Freeze-Out (FO), stages of the reaction are outside the domain of applicability of the fluid dynamical model. After hadronization and FO, the matter is already dilute and can be described well with kinetic models. The initial stages are the most problematic. Non of the theoretical models on the physics market can unambiguously describe the initial stages (see for example  for discussion).

Thus, since we can not unambiguously describe the initial stages of relativistic heavy ion collisions in microscopic models, we have to use some phenomenological models. There are two phenomenological models, which are most frequently discussed in the literature - the Landau and the Bjorken scenarios. These models describe two extreme cases - complete stopping or complete transparency, respectively. Nowadays, experiments have entered the region where the Bjorken model is expected to be applicable, but there is no clear and unambiguous confirmation that experiments follow this scenario. The preliminary experimental results from RHIC do not show transparency - most particle multiplicities (except maybe the most central collisions) do not show a dip in the rapidity spectra, but rather a plateau around mid rapidity , which is a sign of strong stopping. Furthermore, very strong elliptic flow has been measured, which shows a clear peak around mid rapidity . To build such a strong elliptic flow, strong stopping and momentum equilibration are required. Also the $`\overline{p}/p`$ ratio at mid rapidity measured at RHIC  (preliminary) is still far from 1, which tells us that the middle region is not baryonfree.

Our goal is to build a Multi Module Model for ultra-relativistic heavy ion collisions valid for RHIC and LHC energies, and maybe for the most energetic SPS collisions. In the present work only the first step is done - a Two Module Model  \- but a very important one, since the Effective String Rope Model (ESRM)  has been developed for the most problematic module - the module describing the initial stages of collisions.

We describe the initial state with the ESRM; then the hydrodynamical evolution with a QGP EoS is used for the intermediate state. The hydrodymanical calculations are performed with the Los Alamos Particle-in-Cell (PIC) one fluid code . The hydro evolution stops at the FO hypersurface. We present a version of the code assuming that FO (simultaneous chemical and thermal) happens on the simplified toothed hypersurface, where it’s normal vector, $`d\sigma ^\mu `$, is parallel to the flow velocity for every cell. On average this hypersurface approximates the constant time hypersurface. Therefore the flow velocity does not change during the FO process, and the calculations can be done in local rest frame of the matter. Such a surface is also completely timelike; this allows us to avoid the problems discussed in Refs. . The more advanced description of the FO process is planned to be separated into the Third Module.

The EoS presently used in the code was used and presented in Ref. . This phenomenological phase transition model
A) takes the phenomenological EoS for hadronic matter in a simple form, which nevertheless allows one to check different parameterizations discussed in the literature;
B) uses the Bag model EoS to describe QGP;
C) creates a complete EoS, containing pure phases and a region where they coexist, by the Maxwell construction.

The reader may notice that the model is still in a formative stage and a lot of further work remains necessary. We are not yet ready to present the quantitative calculations to be compared with data. Nevertheless we present some preliminary results. First of all, the flow component, $`v_1`$, as a function of rapidity, $`y`$, is presented in Fig. 1, to show that our expectation that the model will generate a third flow component (Refs. ) is reasonable . Our initial state generated by the ESRM indeed produces a strong antiflow in semi-central collisions. The directed $`v_1`$ component appears to be very small, as expected for RHIC energies (for example ). The peaks in $`v_1`$ (actually $`v_3`$) around $`y=\pm 0.3`$ look very high, but this is typical for calculations without thermal smearing . Including thermal smearing will lead to smaller and wider peaks.

We also present in Fig. 2 the evolution of the energy density distribution, $`e(r,t)`$, during the hydrodynamic expansion (notice that $`t=0`$ corresponds to the moment then we start the hydro description, i.e. $`3.5fm/c`$ after the collision). We can see the development of the small but very dense, shock-like domains at the intermediate stages. This effect, at least to some extent, is due to numerical artifacts. The PIC method has a well-known defect called ’ringing’ or marker-clustering . It is called this way because early calculations in two dimensions resulted in density distributions that looked similar to the vibration patterns in sound vibrations - areas of higher and lower densities . Advanced algorithms that mitigate this effect are being implemented in the code.

Another reason which may cause the development of the numerical artifacts is a step-like in beam direction initial energy density distribution (output of ESRM): it has a jump from $`e(x,y)=const`$ inside the matter to $`0`$ in the outside vacuum. In order to avoid (or at least to suppress) the effect of this we are planning to smooth over initial energy density distribution, for example by a Gaussian shape, as proposed in Ref. . Thus it is plausible that the development of the shock-like domains in the simulation is not a physical effect, although we can not rule out other possibilities.

To evaluate the observables a third, FO module is going to be attached to the model. Here the widely used Cooper-Frye model should be essentially improved and modified. First of all conservation laws and the requirement of increasing entropy should be enforced in the module. In addition, particularly for FO across space-like hypersurfaces, realistic, non-equilibrium post FO phase space distributions will have to be used to avoid negative contributions occuring in the naive use of the Cooper-Frye model .


text:
# Quantum fluctuations of coupled dark solitons in trapped Bose-Einstein condensates

## ACKNOWLEDGMENTS

We acknowledge the support from an RGC Earnmarked Grant CUHK4237/01P and the Chinese University of Hong Kong Direct Grant (Grant Nos: 2060148 and 2060150).