text:
# Bulk inhomogeneous phases of anisotropic particles: A fundamental measure functional study of the restricted orientations model

## I Introduction

Onsager first showed that the isotropic-nematic liquid crystal phase transition occurs in systems of anisotropic particles interacting via hard core repulsions Onsager . He studied a system of hard spherocylinders in the limit of infinite anisometry $`\kappa =(L+D)/D\mathrm{}`$ ($`\kappa `$ is the spherocylinder length to breath ratio) using the second virial form of the free energy, which in this limit is exact for the isotropic phase. The effect that higher virial coefficients have in the isotropic-nematic transition was later studied by Zwanzig, who introduced a model of hard prolate uniaxial parallelepipeds with axes oriented along the three perpendicular directions Zwanzig . This peculiar model, which obviously treats the orientational degrees of freedom in an unrealistic way, has the advantage of being accessible to the calculation of higher virial coefficients up to seventh order in the infinite aspect ratio limit. He showed that including all these virials the isotropic-nematic transition also occur, although the exact value of the coexisting nematic density strongly depends on the order of the approximation. The Padé approximant generated by the truncated cluster expansion provides a much more stable sequence of the parameters which characterize the transition Runnels . This stability leaves little room for doubts regarding the existence of the transition in the model. The virial expansion resumed and expressed in the variable $`y=\rho /(1\rho v)`$, with $`\rho `$ the number density of parallelepipeds and $`v`$ their volume, converges more rapidly than the traditional expansion in $`\rho `$, as was shown by Barboy and Gelbart for different hard particle geometries Barboy . Thus, the so-called $`y_3`$ expansion of the Zwanzig model was applied to the study of the isotropic-nematic transition as well as to the study of properties of its interface Moore . For the latter the authors applied the smoothed density approximation of the free energy functional in the spirit of Tarazona’s weighted density approximation for the fluid of hard spheres Tarazona .

The restricted orientations model for hard cylinders was also used to describe the structural properties of molecular fluids near hard walls or confined in a slit. This time the density functional was constructed from the bulk direct correlation function approximated by a linear combination of geometrical functions Rickayzen .

Computer simulations of a variety of models of nonspherical hard core particles showed that the excluded volume effects could not only account for the stability of nematics, but also for the existence of liquid crystal inhomogeneous phases such as the smectic Strobants\_1  and columnar Veerman  phases. Particularly the complete phase diagram of freely rotating hard spherocylinders Bolhuis , including not only smectic, but also a plastic solid phase and different oriented solid phases was calculated. Several density functional theories, all of them based on weighted or modified weighted density approximations, are able to reproduce reasonably well the isotropic-smectic or nematic-smectic transitions Somoza ; Poniewierski ; Velasco  in the whole range of aspect ratios where the smectic is stable, and in some cases, transitions from the isotropic fluid to the plastic or oriented solid phases Graf . In all these approximations the excess free energy is evaluated by integration of the free energy per particle of a reference fluid (typically spheres or hard parallel ellipsoids) evaluated at some weighted or effective density. In some cases, the employed weight is directly the normalized Mayer function between spherocylinders Poniewierski ; Velasco ; in others, it is calculated from the knowledge of the bulk correlation function of the reference fluid Velasco . For the latter case, the term proportional to the Mayer function enters into the integrand as a multiplicative factor of the free energy per particle. The hard sphere free energy functional is recovered in both approaches as the limiting case of $`L=0`$.

The fundamental measure theory (FMT) first developed for hard spheres by Rosenfeld Rosenfeld\_1  was another starting point for constructing a density functional for anisotropic particles. In its general formalism the excess free energy density of the fluid is a function of some weighted densities obtained by convoluting the density profiles with weights which are characteristic functions of the geometry of a single particle whose integrals are the so-called fundamental measures: volume, surface area, and mean radius of the particles. Unfortunately, the Mayer function of two convex anisotropic bodies cannot be decomposed as a finite sum of convolutions of single particle weights Rosenfeld\_2 , which is the keystone for constructing such a functional. Thus, the low density limit of the direct correlation function is no more the Mayer function.

In spite of this, Chamoux and Perera have taken advantage of Rosenfeld’s extension of FMT to hard convex bodies by using it to compute the direct correlation function and patching out the low density limit with the exact Mayer function Chamoux\_1 . In this way they have obtained the equation of state for various convex hard bodies (such as hard ellipsoids, spherocylinders, and cut spheres), have predicted ordered phases and, recently, have study demixing in binary mixtures of rigidly ordered particles Chamoux\_1 ; Chamoux\_2 .

Following a similar procedure a density functional for anisotropic particles has been proposed which interpolates between the Rosenfeld’s hard sphere functional and Onsager’s functional for elongated rods. The resulting model was tested by calculating the isotropic-nematic transition in systems of hard spherocylinders and hard ellipsoids Cinacchi .

Although the authors of this work suggest that the resulting theory can be applied to the study of inhomogeneous systems, the huge computational efforts that their numerical implementations involve is the reason for the absence of any result in this direction. One way to circumvent this difficulty is to reduce the continuous orientational degree of freedom to three discrete orientations (Zwanzig model). Implementing this idea some authors have recently applied the Zwanzig model to the study of interfacial properties of the hard rod fluid interacting with a hard wall or confined in a slit, for a one-component vanRoij  and a polydisperse mixture Martinez-Raton\_1 , and also to the study of bulk and interfacial properties of hard platelet binary mixtures Harnau . All these models are based on Onsager’s density functional approximation. The increase of the number of allowed orientations in this functional particularized for hard spherocylinders results in the presence of an artificial nematic-nematic transition in the one component fluid as the authors of Ref. vanRoij\_1  have shown. This result indicates that certain cares must be taken in the direct extrapolation of the results obtained from this theory.

FMF was also constructed for a mixture of parallel hard cubes combining Rosenfeld’s original ideas with a dimensional cross over constraint Cuesta . The latter appears to be very important to describe correctly the structure of inhomogeneous fluids in situations of high confinement and to describe well the structural properties of the solid phase Tarazona\_1 . The dimensional cross over has been used as an important ingredient to develop a density functional for a binary mixture of hard spheres and needles, assuming that the needles are too thin to interact with each other directly Schmidt .

Taking a ternary mixture of parallel hard cubes and scaling each species along one of the three Cartesian axes with the same scaling factor a FMF for the Zwanzig model is obtained. This functional has already been applied to the study of the effect that polydispersity has on the stability of the biaxial phase in a binary mixture of rods and plates Martinez-Raton\_2  and on the relative stability of the smectic and columnar phases due to the presence of polydispersity Martinez-Raton\_3 .

The FMF for Zwanzig’s model in the homogeneous limit coincide with the scaled particle theory and thus with the so-called $`y_3`$ expansion which, as pointed out before, first began to be used in Ref. Moore  as a model to study the isotropic-nematic phase transitions in fluids of hard parallelepipeds. But this functional, through its minimization, also allows us to calculate inhomogeneous density profiles. This functional has been applied recently to study the isotropic-nematic interface of a binary mixture of hard platelets Bier . Its structural and thermodynamic properties resulting from the FMF minimization show complete wetting by a second nematic. The same phenomenon was found in a binary mixture of hard spherocylinders Shundyak .

The phase diagram for Zwanzig’s model including the smectic, columnar, and solid phases has never been carried out, only spinodal instability boundaries have been traced Martinez-Raton\_3 . The main purpose of this work is to obtain the complete phase diagram for this model and to compare the results with the only existing simulation of the lattice version of the model, which has been carried out for two different aspect ratios Casey . This will test the predictive power of the FMF for anisotropic inhomogeneous phases. As a particular case, the system of parallel hard cubes will be studied. In Ref. Martinez-Raton\_4  a bifurcation analysis and a Gaussian parametrization of the density profiles were used to calculate the free energy and pressure of the solid phase. Here a free minimization will be performed to calculate not only the solid but also the smectic and columnar phases and compare the obtained results with recent simulations of parallel hard cubes Jagla ; Groh .

## II FMF for Zwanzig model

The FMF for hard parallelepipeds was already described in detail in Ref. Cuesta . A brief summary of the theory will be presented here putting emphasis on its numerical implementation to calculate the equilibrium inhomogeneous phases.

A ternary mixture of hard parallelepipeds of cross section $`\sigma ^2`$ and length $`L`$ with their uniaxial axes pointing to the $`x`$, $`y`$, or $`z`$ directions is described in terms of their density profiles $`\rho _\mu (𝐫)`$ ($`\mu =x,y,z`$). Following the FMT for hard parallelepipeds in three dimensions the excess free energy density in reduced units can be written as Cuesta

$`\mathrm{\Phi }_{\mathrm{exc}}(𝐫)=\mathrm{\Phi }^{(1)}(𝐫)+\mathrm{\Phi }^{(2)}(𝐫)+\mathrm{\Phi }^{(3)}(𝐫),`$ (1)
where the $`\mathrm{\Phi }^{(k)}`$’s are

$`\mathrm{\Phi }^{(1)}`$ $`=`$ $`n_0\mathrm{ln}(1n_3),`$ (2)
$`\mathrm{\Phi }^{(2)}`$ $`=`$ $`{\displaystyle \frac{𝐧_1𝐧_\mathrm{𝟐}}{1n_3}},`$ (3)
$`\mathrm{\Phi }^{(3)}`$ $`=`$ $`{\displaystyle \frac{n_{2x}n_{2y}n_{2z}}{(1n_3)^2}},`$ (4)
with weighted densities

$`n_\alpha (𝐫)={\displaystyle \underset{\mu }{}}\left[\rho _\mu \omega _\mu ^{(\alpha )}\right](𝐫),`$ (5)
i.e., they are sums of convolutions of the density profiles with the following weights:

$`\omega _\mu ^{(0)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{8}}{\displaystyle \underset{k=1}{\overset{3}{}}}\delta \left({\displaystyle \frac{\sigma _\mu ^k}{2}}|x_k|\right),`$ (6)
$`\omega _\mu ^{(3)}(𝐫)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{3}{}}}\theta \left({\displaystyle \frac{\sigma _\mu ^k}{2}}|x_k|\right),`$ (7)
$`\omega _\mu ^{(1j)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{2\theta \left(\frac{\sigma _\mu ^j}{2}|x_j|\right)}{\delta \left(\frac{\sigma _\mu ^j}{2}|x_j|\right)}}\omega _\mu ^{(0)}(𝐫),`$ (8)
$`\omega _\mu ^{(2j)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{\delta \left(\frac{\sigma _\mu ^j}{2}|x_j|\right)}{2\theta \left(\frac{\sigma _\mu ^j}{2}|x_j|\right)}}\omega _\mu ^{(3)}(𝐫),`$ (9)
where the notation $`x_k`$ ($`k=1,2,3`$) for the $`x`$, $`y`$, and $`z`$ coordinates has been employed. The functions $`\delta (x)`$ and $`\theta (x)`$ are the usual delta Dirac and Hevisaide functions and $`\sigma _\mu ^j=\sigma +(L\sigma )\delta _\mu ^j`$ with $`\delta _\mu ^j`$ the Kronecker delta.

The following constraints on the density profiles were imposed: (i) The solid phase has the simple parallelepipedic unit cell with uniaxial symmetry, i.e., the periods in the three spatial directions are $`d_{}`$ for $`x,y`$ and $`d_{}`$ for $`z`$. The orientational director is selected parallel to $`z`$. (ii) The density profile of each species has the form

$`\rho _\mu (𝐫)=\rho \gamma _\mu {\displaystyle \underset{𝐤=\mathrm{𝟎}}{\overset{𝐧}{}}}\alpha _𝐤^{(\mu )}{\displaystyle \underset{j=1}{\overset{3}{}}}\mathrm{cos}\left(q_jk_jx_j\right),`$ (10)
where $`q_j=2\pi /d_j`$ is the wave number along the $`j`$ direction, $`𝐤=(k_1,k_2,k_3)`$ is the vector defined by the reciprocal lattice numbers, and $`𝐧=(n_1,n_2,n_3)`$ is the vector at which the harmonic expansion is truncated. Thus, Eq. (10) is the Fourier expansion of the density profiles $`\rho _\mu (𝐫)`$ truncated at some $`𝐧`$. This cutoff is selected in such a way that it guarantees small enough values of $`\alpha _𝐧^{(\mu )}`$. The first Fourier amplitudes of all species are fixed to one ($`\alpha _\mathrm{𝟎}^{(\mu )}=1`$) and consequently $`V_{\mathrm{cell}}^1_{V_{\mathrm{cell}}}𝑑𝐫\rho _\mu (𝐫)=\rho \gamma _\mu `$ with $`V_{\mathrm{cell}}=d_{}^2d_{}`$ the unit cell volume, $`\rho `$ the mean total density over the unit cell, and $`\gamma _\mu `$ the occupancy probability of species $`\mu `$ in the unit cell, which obviously fulfills the condition $`_\mu \gamma _\mu =1`$.

In the plastic solid phase these occupancy probabilities are $`1/3`$ for each species while they deviate from this value in the oriented solid phase. The uniaxial symmetry also implies that $`\gamma _x=\gamma _y=\gamma _{}`$, $`\gamma _z=\gamma _{}=12\gamma _{}`$ and $`\rho _x(x,y,z)=\rho _y(y,x,z)`$, $`\rho _z(x,y,z)=\rho _z(y,x,z)`$. Thus, the Fourier amplitudes verify $`\alpha _{(k_1,k_2,k_3)}^{(x)}=\alpha _{(k_2,k_1,k_3)}^{(y)}`$ and $`\alpha _{(k_1,k_2,k_3)}^{(z)}=\alpha _{(k_2,k_1,k_3)}^{(z)}`$. The total number of Fourier amplitudes \[except the $`(0,0,0)`$ term of all species\] is reduced by these symmetries to $`N_\alpha =(n_{}+1)(n_{}+1)(3n_{}+4)/22`$, ($`n_1=n_2n_{}`$, $`n_3n_{}`$) independent variables. These variables together with $`\gamma _{}`$, $`q_{}`$ and $`q_{}`$ span the variable space in which the FMF must be minimized.

The density profiles of columnar and smectic phases are obtained from Eq. (10) substituting $`𝐧=(n_{},n_{},0)`$ and $`𝐧=(0,0,n_{})`$. From the definition (5), Eqs. (6)-(9) and the density profiles (10), the weighted densities can be easily calculated resulting in

$`n_\alpha (𝐫)`$ $`=`$ $`\rho {\displaystyle \underset{\mu ,𝐤}{}}\gamma _\mu \alpha _𝐤^{(\mu )}\chi _{\alpha ,𝐤}^{(\mu )}{\displaystyle \underset{j=1}{\overset{3}{}}}\mathrm{cos}\left(q_jk_jx_j\right),`$ (11)
$`\chi _{0,𝐤}^{(\mu )}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{3}{}}}\varphi _0\left(\xi _{j,𝐤}^{(\mu )}\right),`$ (12)
$`\chi _{3,𝐤}^{(\mu )}`$ $`=`$ $`v{\displaystyle \underset{j=1}{\overset{3}{}}}\varphi _3\left(\xi _{j,𝐤}^{(\mu )}\right),`$ (13)
$`\chi _{1j,𝐤}^{(\mu )}`$ $`=`$ $`\sigma _j^\mu {\displaystyle \frac{\varphi _3\left(\xi _{j,𝐤}^{(\mu )}\right)}{\varphi _0\left(\xi _{j,𝐤}^{(\mu )}\right)}}\chi _{0,𝐤}^{(\mu )},`$ (14)
$`\chi _{2j,𝐤}^{(\mu )}`$ $`=`$ $`{\displaystyle \frac{1}{\sigma _j^\mu }}{\displaystyle \frac{\varphi _0\left(\xi _{j,𝐤}^{(\mu )}\right)}{\varphi _3\left(\xi _{j,𝐤}^{(\mu )}\right)}}\chi _{3,𝐤}^{(\mu )},`$ (15)
with $`v=L\sigma ^2`$ the particle volume, $`\varphi _0(x)=\mathrm{cos}x`$, $`\varphi _3(x)=\mathrm{sin}x/x`$, and $`\xi _{j,𝐤}^{(\mu )}=q_jk_j\sigma _j^\mu /2`$.

The substitution of Eqs. (10) and (11) into the free energy per unit cell

$`\mathrm{\Phi }{\displaystyle \frac{\beta }{V_{\mathrm{cell}}}}`$ $`=`$ $`V_{\mathrm{cell}}^1{\displaystyle _{V_{\mathrm{cell}}}}𝑑𝐫\left[\mathrm{\Phi }_{\mathrm{id}}(𝐫)+\mathrm{\Phi }_{\mathrm{exc}}(𝐫)\right],`$ (16)
$`\mathrm{\Phi }_{\mathrm{id}}(𝐫)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\rho _\mu (𝐫)\left[\mathrm{ln}\left(\rho _\mu (𝐫)\mathrm{\Lambda }_\mu ^3\right)1\right],`$ (17)
with $`\mathrm{\Phi }_{\mathrm{id}}(𝐫)`$ the ideal part of the free energy density, and its minimization with respect to the $`N_\alpha +3`$ variables allows the calculation of the equilibrium free energy and the density profiles of inhomogeneous phases.

To characterize the structure and orientational order of these phases the following total density and order parameter profiles will be used:

$`\rho (𝐫)`$ $`=`$ $`{\displaystyle \underset{\mu }{}}\rho _\mu (𝐫),`$ (18)
$`Q(𝐫)`$ $`=`$ $`1{\displaystyle \frac{3}{2}}{\displaystyle \frac{\left[\rho _x(𝐫)+\rho _y(𝐫)\right]}{\rho (𝐫)}}.`$ (19)
The selection of $`Q(𝐫)`$ as an order parameter is motivated by its uniaxial symmetry property $`Q(x,y,z)=Q(y,x,z)`$ and its uniform limit value $`Q=13\gamma _{}`$ ($`1/2Q1`$), which coincides with the usual definition of the nematic order parameter for the Zwanzig model: $`Q=0`$ ($`\gamma _{}=1/3`$) for the isotropic phase and $`Q=1`$ ($`\gamma _{}=0`$) for the perfectly aligned nematic phase. Although the solid and columnar phases might have local biaxiality \[$`\rho _x(x,y,z)\rho _y(x,y,z)`$\], the integral over the unit cell of any previously defined biaxial order parameter is always equal to zero as a consequence of the symmetries of the density profiles.

## III Phase diagrams

The phase diagrams presented in this work were calculated for a set of aspect ratios ranging from $`\kappa =0.1`$ to $`\kappa =10`$, corresponding to the aspect ratios of the most anisotropic oblate and prolate parallelepipeds studied here. The volume of all particles (cubes or prolate or oblate parallelepipeds) are fixed to 1 and thus the mean packing fraction $`\eta `$ is equal to the mean density $`\rho `$. From the equation $`v=L\sigma ^2=1`$ the parallelepiped edge lengths $`L`$ and $`\sigma `$ can be calculated as a function of the aspect ratio $`\kappa =L/\sigma `$ as $`L=\kappa ^{2/3}`$ and $`\sigma =\kappa ^{1/3}`$. For each $`\kappa `$, fixing the mean density $`\rho `$ and using appropriate initial guesses for the $`N_\alpha +3`$ variables with symmetries corresponding to the smectic, columnar or solid phases, the energy per unit cell (16) was minimized and thus the free energy for each phase was obtained. Varying $`\rho `$ and repeating the former steps the free-energy branches of the different inhomogeneous phases have been calculated. The common tangent construction allowed the calculation of the coexisting densities between those phases in the case of first order transitions. To evaluate numerically the three dimensional integral of the free energy density (16) a Gauss-Chebyshev quadratures has been employed.

### III.1 Parallel hard cubes

This subsection is devoted to the study of the parallel hard cube system ($`\kappa =1`$). The PHC equation of states of the fluid and solid phases as obtained from the FMT and the Monte Carlo simulation results are compared. While the solid phase is very well described with this formalism the exact location of the fluid-solid transition is very poorly estimated. The fundamental reasons of this difference are discussed here through a critical analysis of the fluid equation of state resulting from the FMF. It will be shown that possible modifications of the FMF slightly improve the location of the transition point at the expense of the correct description of the solid branch.

In Ref. Martinez-Raton\_2  the PHC fluid was already studied with the same FMF but using a Gaussian parametrization for the density profile. Through a minimization procedure and also from a bifurcation analysis a second-order fluid-solid transition was found at $`\rho =0.3143`$ with a lattice period $`d=1.3015`$ and a fraction of vacancies $`\nu =0.3071`$ Martinez-Raton\_2 . Recent simulations on the same system also showed a second-order transition to the solid but with very different transition parameters $`\rho =0.48\pm 0.02`$ in Ref. Jagla  and $`\rho =0.533\pm 0.010`$ in Ref. Groh . No evidence for the vacancies predicted by FMT was found, although the authors recognized that the vacancies might be suppressed by the boundary conditions in the small systems accessible to simulations Groh .

The main problem of the FMF for hard cubes is that it recovers in the homogeneous limit the scaled-particle equation of state, which overestimates the pressure calculated from the exact virial expansion up to seventh order. This expansion has a maximum at $`\rho 0.6`$ and then goes down very quickly to reach negative values Swol . The poorly convergent character of the virial series makes it impossible to construct an equation of state for hard cubes, such as the Carnahan-Starling equation for the hard-sphere fluid, which estimates reasonably well all the known virial coefficients and diverges at close packing. On the other hand, it is well known that the FMF describes accurately the fluid structure in situations of high confinement, including the solid phase near close packing. For example, at high densities the functional recovers the cell theory, which is asymptotically exact when the packing fraction goes to 1, and also compares reasonably well with computer simulations Groh . These nice properties are a consequence of a fundamental restriction, namely, the dimensional cross-over Cuesta , imposed in the construction of the FMF. The latter implies that the functional in dimension $`D`$ reduces to the functional in dimension $`D1`$ when the original density profile is constrained to $`D1`$ dimensions, i.e., $`\rho ^{(D)}(𝐫)=\rho ^{(D1)}(𝐫)\delta (x_D)`$, where $`x_D`$ is the coordinate that is eliminated on going from $`D`$ to $`D1`$ dimensions.

One possible procedure to improve the description of the uniform fluid of hard cubes at the level of the FMF is to follow the same method used in Refs. Tarazona\_2  and Roth , in which the hard-sphere Carnahan-Starling equation of state is imposed through the modification of the third term $`\mathrm{\Phi }^{(3)}`$ \[see Eq. (1)\] of the excess free-energy density while keeping the exact density expansion of the direct correlation function up to first order. Unfortunately the absence of a good equation of state for the PHC fluid with the already mentioned properties makes this procedure less systematic compared to that of hard spheres Tarazona\_2 ; Roth .

Following this purpose the original excess free-energy density for hard cubes (1) is now substituted by

$`\mathrm{\Phi }_{\mathrm{exc}}(𝐫)={\displaystyle \underset{k=1}{\overset{3}{}}}f_k\left(n_3(𝐫)\right)\mathrm{\Phi }^{(k)}(𝐫),`$ (20)
with the $`f_k(n_3)`$’s selected in such a way as to keep the correct first order density expansion of the direct correlation function and to obtain the right virial expansion up to the seventh order of the PHC equation of state. As the original FMF for hard cubes gives the third virial coefficient correctly, these conditions imply that $`f_{1,2}(n_3)1+O(n_3^2)`$ and $`f_3(n_3)1+O(n_3)`$ for small $`n_3`$. Two further important conditions imposed on the $`f_k(n_3)`$’s are their limiting behavior when the local packing fraction tends to unity: $`lim_{n_31}f_\alpha (n_3)=1`$, which asymptotically guarantees the correct cell-theory limit, and the positive signs of their values, which guarantee the convexity of the fluid free energy. Unfortunately this procedure breaks the dimensional cross-over property, but in principle should describe the fluid-solid transition in hard cubes better.

Among all the functions $`f_\alpha `$’s that have been tried, even those which give better results \[the particular case of $`f_{1,2}(n_3)=1`$\] are far from getting the transition point near the simulation one. In Fig. 1(a) the scaled-particle equation of state, the improved equation of state

$`\beta P=\rho +\rho {\displaystyle \frac{\mathrm{\Phi }_{\mathrm{exc}}}{\rho }}\mathrm{\Phi }_{\mathrm{exc}},`$ (21)
with $`\mathrm{\Phi }_{\mathrm{exc}}`$ being the uniform limit of Eq. (20), and finally, the symmetric Padé approximant of the seventh-order virial series are plotted.

In the first two curves the bifurcation points are shown. The new bifurcation point calculated from Eq. (20) is located at $`\rho =0.3378`$, and the period and fraction of vacancies of the solid are $`d=1.3249`$ and $`\nu =0.2143`$. As can be seen from Fig. 1(a), the new equation of state still overestimates the fluid pressure, but to a lesser extent. Although the new functional gets a higher transition density and the fraction of vacancies decreases, there is still disagreement between theory and simulations. The equation of state of the PHC solid calculated from the minimization of the original FMF with respect to the Gaussian density profiles compare very well with simulations for densities $`\rho 0.5`$, whereas the modified version underestimates the solid pressure.

At this point the main conclusion that can be drawn is that the modification of the FMF in order to improve the description of the uniform fluid spoils the good description of the solid phase. As the modification of the FMF was done at the expense of loosing the dimensional cross-over property (and this spoils the good description of highly inhomogeneous phases), and the modified versions do not show too many differences in the prediction of the fluid-solid transition, it is worthless to use them to study nonuniform phases.

Setting $`q_{}=q_{}=q=2\pi /d`$, $`\gamma _\mu =1/3`$, and $`\alpha _𝐤^{(\mu )}=\beta _𝐤`$ in the density profiles (10) and minimizing the FMF, Eq. (16), of parallel hard cubes ($`\kappa =1`$) with respect to the Fourier amplitudes and the wave number $`q`$, the free energy per unit cell for solid \[$`𝐧=n(1,1,1)`$\], columnar \[$`𝐧=n(1,1,0)`$\] and smectic \[$`𝐧=n(0,0,1)`$\] phases were obtained. The results are shown in Fig. 2. From the isotropic liquid at the same density $`\rho =0.3134`$ three inhomogeneous solutions: solid, columnar, and smectic, bifurcate, with the solid phase being the stable one. While the free energy difference between solid and columnar phases is relatively small, the smectic phase is clearly thermodynamically unfavorable.

The number of Fourier amplitudes necessary to describe adequately the density profile increases with the density, and thus the numerical calculations becomes more and more time consuming. Nevertheless, the scenario shown in Fig. 2, with the solid being the only stable phase, occurs at high densities as the simulations and cell-theory have confirmed Groh . The minimization of the FMF using a Gaussian parametrization of the density profiles of columnar and solid phases shows very similar quantitative results Groh . In fact the equation of state of the parallel hard-cube solid from FMT calculations with this parametrization compares very well at high densities with simulations Martinez-Raton\_3 . The results presented here are much more accurate than those obtained through the Gaussian parametrization.

### III.2 Prolate parallelepipeds

This subsection is devoted to study the phase diagram of prolate particles ($`\kappa >1`$). The results obtained from numerical minimization of the FMF of parallelepipeds with fixed $`\kappa =4.5`$ are shown in Fig. 3. The free energies per unit volume of those phases which are stable in some range of densities are plotted. As can be seen the isotropic phase undergoes a first-order phase transition to the so-called discotic smectic (DSm) phase. This peculiar phase is a layered phase (similar to the smectic phase) but with the long axes of the parallelepipeds lying within the layers. There is no orientational order in the layers, what means that the order parameter $`Q(z)`$ reaches negative values at the positions of the density peaks. The density and order parameter profiles of the DSm phase at $`\rho =0.3`$ are plotted in Fig. 4.

The period in units of the small particle length is $`d/\sigma =1.2796`$ which means that the particles with long axes perpendicular to the layers (preferentially localized at the center of the interlayer space) intersect approximately three adjacent layers.

Simulations of the Zwanzig model with $`\kappa =5`$ on a lattice showed an I-DSm transition at a density between 0.47 and 0.55 Casey . Although the results were obtained for a lattice spacing of $`1/3`$ (in units of the shortest particle dimension) the simulations were repeated for values $`1/9`$ and $`1/27`$ without changes in the stability of the DSm phase. Thus, the authors concluded that this layered phase may persist in the continuum limit Casey . The difference in the transition density found from FMT ($`0.2868`$) and from simulations ($`0.5`$) can be explained using two arguments: (i) As was already pointed out in Sec. II, the FMF in the uniform density limit considerably overestimates the isotropic fluid pressure and thus the theory underestimate the transition densities between homogeneous and inhomogeneous phases. (ii) The transition densities should decrease upon decreasing the lattice spacing in simulations, as the results for the freezing of parallel hard cubes on a lattice (occurring at $`\rho =0.568`$ for an edge length equal to two lattice spacings, at $`\rho =0.402`$, for six lattice spacings, and at $`\rho =0.314`$ for the continuum) illustrate Lafuente .

Increasing the mean density further the DSm phase undergoes a first-order transition to the columnar phase as Fig. 3 shows. The restriction of parallelepiped orientations enhances the columnar phase stability even for prolate particles as a phase diagram, to be described below, will show. This phenomenon can be understood if the Zwanzig model is interpreted as a ternary mixture of particles. Simulations on a binary mixture of parallel spherocylinders with different aspect ratios (specifically 2 and 2.9) show that, instead of a continuous nematic-smectic transition typical of the pure component system, the mixture exhibits a first-order nematic-columnar phase transition Stroobants\_2 . This result was explained by the poorer packing of rods of different lengths in the smectic phase as compared to that of rods of the same length. Simulations and theory show that one of the most important effects that the aspect ratio polydispersity has on the phase behavior of hard spherocylinders Polson  and binary mixtures of oblate and prolate particles vanKoij ; Martinez-Raton\_3  is the enhancement of the columnar phase stability. All these results show that the columnar phase can be stable even for mixture of particles with different shapes. Although the constituent particles of the Zwanzig model have the same shape, the restriction of their orientations changes strongly its relative packing and thus for some $`\kappa `$’s enhance the columnar phase stability with respect to other phases.

At higher density the columnar phase exhibits a continuous phase transition to an oriented solid phase of prolate parallelepipeds, as shown in Fig. 5(a). The density and order parameter profiles of the columnar phase at the bifurcation point ($`\rho =0.3748`$) are shown in Fig. 6. The periods of the solid phase along the perpendicular and parallel directions are $`d_{}/\sigma =1.2690`$ and $`d_{}/L=1.5170`$, respectively. From the equation $`\rho =(1\nu )V_{\mathrm{cell}}^1`$ ($`V_{\mathrm{cell}}=d_{}^2d_{}`$ being the unit cell volume), the fraction of vacancies of the solid $`\nu `$ can be calculated as $`0.0845`$. The continuous nature of the columnar-oriented solid transition changes to first order at some $`\kappa `$ between $`4`$ and $`4.5`$, as Fig. 5(b) shows for $`\kappa =4`$.

The order parameter $`Q(𝐫)`$ is very high in the unit cell except in its borders, where it exhibits small oscillations \[see Fig. 6(b)\]. These oscillations are a consequence of the microsegregation of species “$`x`$” and “$`y`$” in the newly formed solid phase, which is preferentially formed by particles of species “$`z`$” localized around the position $`(x^{},y^{})=(0,0)`$.

This feature is shown in Fig. 7, where the sum of the density profiles of species “$`x`$” and “$`y`$” \[$`\rho _{}(𝐫)=\rho _x(𝐫)+\rho _y(𝐫)`$\] is plotted. While the columnar packing is responsible for the presence of the local maxima at the center of the unit cell, the species “$`x`$” and “$`y`$” begin to segregate to the borders of the cell $`(\pm 0.5,0)`$ and $`(0,\pm 0.5)`$, respectively (see the four local maxima at these positions) as the new solid phase is formed. The long axes of the perpendicular species lie on the lateral sides of the parallelepipedic unit cell, while their centers of mass are preferentially localized at the centers of these sides.

The calculation of the free-energy branches for several stable inhomogeneous phases and the phase transitions between them (as it was described for $`\kappa =4.5`$) has been carried out for 15 values of $`\kappa `$ (ten of them in the range $`1\kappa 5`$ and five of them in the range $`5\kappa 10`$). The resulting phase diagram is plotted in Fig. 8. The isotropic phase of prolate parallelepipeds with $`1\kappa 3.5`$ undergoes a continuous phase transition to the plastic solid phase. The transition points are joined with the spinodal line that has been calculated through the divergence of the structure factor. Notwithstanding that a functional minimization was carried out for each $`\kappa `$ to check the continuous nature of the transitions. The plastic solid is stable for $`\kappa 2.5`$ up to densities around $`0.5`$. At these values the numerical minimization turns out to be cumbersome because of the large number of Fourier amplitudes necessary to correctly describe the inhomogeneous profiles. Thus, the high density part of the phase diagram ($`\rho 0.5`$) has not been calculated with the numerical procedure described above. At higher densities a Gaussian-type parametrization of the density profiles is required, which obviously has a lower degree of accuracy.

For $`\kappa =2.95`$ the plastic solid exhibits a very weak first-order phase transition to the discotic smectic phase (labeled as 1 in Fig. 8), and the latter a phase transition to the columnar phase at higher densities. But the most representative region of the phase diagram where the discotic smectic is stable is for $`\kappa `$ around 4.5 where this layered phase exhibits a first-order phase transition to columnar phase (the shaded area of Fig. 8 limits the instability region against phase separation between both phases). For $`\kappa `$ between 4 and 5 the columnar phase undergoes a phase transition (first order for $`\kappa =4`$ and continuous for other values shown) to the oriented solid phase.

The nematic phase begins to be stable for $`\kappa >5`$ with its stability region bounded below by the first order isotropic-nematic transition and above by a continuous nematic-smectic transition. Finally, the smectic region is bounded above by a continuous transition to the oriented solid phase (see Fig. 8).

Again the nematic-smectic transition points are joined with spinodal lines and for each $`\kappa `$ a minimization was carried out to check numerically the continuous character of the transition (the smectic solution begins to be stable right at the spinodal). In Ref. Martinez-Raton\_5  a bifurcation analysis with the same functional was carried out to study the nature of the nematic-smectic transition. A thermodynamic and mechanical stability analysis showed that the nematic-smectic transition is first order, which is in contradiction with the numerical minimization results presented here. A possible reason that justifies this contradiction could be that the $`N`$-Sm transition is very weakly first order, so weak that the numerical accuracy used in the functional minimization can not decide about its nature. Another possibility is that the numerical accuracy failure is somewhere in the bifurcation analysis. A careful revision of this analysis is certainly called for in order to settle this point.

The available simulation results for freely rotating hard spherocylinders show that the isotropic phase exhibits a transition to the solid phase for $`0\kappa 4.1`$ (the solid is plastic for $`\kappa 1.35`$ and oriented for $`1.35\kappa 4.1`$) while the isotropic-smectic and nematic-smectic transitions begin at $`\kappa =4.1`$ and $`4.7`$, respectively Bolhuis  \[notice that for hard spherocylinders the length-to-breadth ratio is $`\kappa =(L+D)/D`$\]. We can see that, despite the different particle geometry and the restricted orientations of the Zwanzig model, the agreement for the threshold $`\kappa `$ at which spatial instabilities to the solid and smectic phase destabilize the homogeneous phases is rather good. Also the qualitative picture is similar: elongated rods form smectics, and more symmetric particles form solids. The main difference between them is that the Zwanzig phase diagram presents regions where the columnar and discotic smectic phases are stable, a difference due to the restriction of orientations.

### III.3 Oblate parallelepipeds

The phase diagram of oblate parallelepipeds ($`\kappa <1`$) is shown in Fig. 9. The main differences after comparing the phase diagrams of prolate (Fig. 8) and oblate particles are that in the latter: (i) The smectic is no more a stable phase. (ii) The region of columnar phase stability is considerably larger. (iii) The stability region of the plastic solid is reduced (in fact this phase is stable only up to $`\kappa ^12.5`$) at the expense of that of the discotic smectic phase. (iv) The transitions to the latter are strongly first order in nature (except for $`\kappa ^1=4.5`$). (v) The oriented solid phase is replaced by a perfectly oriented solid in which “$`x`$” and “$`y`$” species are absent. This phase, after scaling in the $`z`$ direction, is the same as the solid of parallel hard cubes. A solution from the FMF minimization with three dimensional spatial modulations and with $`\gamma _{}0`$ has not been found in the parallelepipedic unit cell (the case of face-centered or body-centered cubic unit cells have not been tried here).

Finally, similar by to what happens with prolate parallelepipeds, the nematic phase begins to be stable at $`\kappa ^15`$. It undergoes a continuous phase transition to the columnar phase (the transition points of Fig. 9 are joined with the spinodal curve).

The parallelepipeds with $`\kappa ^1=1.5`$ exhibit an interesting phase behavior. At low densities the isotropic phase destabilizes with respect to the columnar phase and not with respect to the PS phase. This example shows that the prediction for phase transitions using only the spinodal instability calculations can generate uncertainties about the possible symmetries of the inhomogeneous phases. In fact these calculations do not allow to decide in this example if the new phase is a plastic solid or a columnar phase. Only by a complete minimization of the FMF, could it be concluded that the columnar phase is the stable one.

In Fig. 10 (a) the density profiles of perpendicular \[$`\rho _{}(z)`$\] and parallel \[$`\rho _{}(z)`$\] species are shown for the discotic smectic phase of oblate particles with $`\kappa ^1=2.5`$, while the order parameter profile is shown in Fig. 10(b). This discotic smectic phase coexists at $`\rho =0.4244`$ with the plastic solid phase. The period in units of the large parallelepiped edge length is $`d/\sigma =1.2142`$. The random orientation of the uniaxial axes within the layers is confirmed by the high negative values of the order parameter at the position of the density peak of the perpendicular species. The main difference between the DSm of Fig. 10 and that of Fig. 7 is that the “z” species is now localized preferentially not at the center of the interlayer space, but near the smectic layers \[see Fig. 10(a)\], exhibiting two local maxima at each side of the layer. This effect can be explained by the depletion force that the perpendicular species exerts on the parallel one.

The $`N`$-Sm ($`NC`$) and the Sm-OS ($`C`$-POS) transition lines of Figs. 8 and 9 converge asymptotically to $`\rho =0.3143`$, the value of the fluid-solid bifurcation density, as $`\kappa \mathrm{}`$ ($`\kappa ^1\mathrm{}`$). The reason for this is that upon increasing $`\kappa `$ ($`\kappa ^1`$) the number of rods (plates) with orientation perpendicular to the director becomes vanishing small, and then the system is, after rescaling the $`z`$ direction, almost equivalent to a system of parallel cubes.

Simulations of the Zwanzig model on a $`15\times 15\times 15`$ lattice with spacing $`1/3`$ show that oblate parallelepipeds with dimensions $`5\times 5\times 1`$ undergo a transition to a phase exhibiting a columnar structure Casey  at a density somewhere between $`0.55<\rho <0.65`$. On increasing the system size to $`30\times 30\times 30`$ the global columnar order disappears, but local correlations persist in the fluid with particle alignment distributed evenly among the three available orientations. In the same work the effect that orientational constraints have on the stability of the inhomogeneous phases was studied. While a system of biaxial $`5\times 3\times 1`$ “tiles” without orientational constraints (except those inherent to the Zwanzig model) stabilizes in a smectic phase with the shortest edge lengths perpendicular to the layers, the system composed by “tiles” with all their long edge lengths parallel to each other exhibits a phase transition to the smectic phase with these edge lengths perpendicular to the layers, similar to what is found here for uniaxial oblate parallelepipeds (the discotic smectic phase).

Simulations of hard cut spheres show that for $`\kappa =0.3`$ there is an isotropic-solid transition, for $`\kappa =0.2`$ an isotropic-columnar transition (the isotropic phase might instead be a peculiar “cubatic” phase) and for $`\kappa =0.1`$ a nematic-columnar one Veerman . From these results it can be concluded that the effect that the degree of particle anisotropy has on the symmetry of the stable phases for both cut spheres and hard parallelepipeds with restricted orientations, is qualitatively similar.

## IV Conclusions

The goal of this article has been the calculation of the phase diagram of the Zwanzig model for prolate and oblate parallelepipeds centering the attention on the phase transitions to inhomogeneous phases. For this purpose the fundamental measure functional for hard parallelepipeds with restricted orientations has been used. This functional is exactly the same for any particle shape (prolate and oblate depending on $`\kappa `$), which allows for a unified study of the phase behavior of both kinds of particles. A free minimization of the functional was carried out with the only constraints of choosing a parallelepipedic unit cell and of imposing uniaxial symmetry in the inhomogeneous phases. The latter is justified by uniaxial symmetry of the particles. The degree of approximation to the exact density profiles was controlled by the cutoff imposed on the reciprocal lattice numbers in the Fourier expansion of the density profiles.

A system of parallel hard cubes was separately studied, which was motivated by recent simulations on this system Jagla ; Groh . Applying a modified versions of the FMF to improve the description of the PHC liquid, along the same lines as Refs. Tarazona\_2  and Roth , the continuous transition point to the solid phase and the equation of state of the solid were calculated from the divergence of the structure factor and from the functional minimization with respect to Gaussian density profiles. Although the transition density and fraction of vacancies change in the right direction, these results are still far from the simulations. In fact, the solid phase is poorly described by the new functional. The poor convergence of the PHC virial series does not make this procedure as effective as for hard spheres. Further refinement of the method and the proper inclusion of vacancies in simulations of the solid phase will probably improve the agreement between theory and simulations.

The original FMF for PHC was minimized to study the relative stability of the smectic, columnar, and solid phases, starting at low densities from the bifurcation point. The solid phase is the only stable phase, followed by the columnar and the smectic (in order of energy stability). At high densities the same behavior is shown from calculations using cell theory, functional minimization with Gaussian density profiles and computer simulations Groh .

The system of prolate and oblate parallelepipeds exhibits a very rich phase behavior. Apart from the plastic or oriented solid, smectic, and columnar phases, which are present also in systems of prolate (spherocylinders Bolhuis ) and oblate (cut spheres Veerman ) particles, a new phase appears: the discotic smectic, the existence of which was confirmed by simulations Casey . The close relation between the particle anisotropy and symmetry of the stable phases (elongated particles form smectics, flattened one form columnars and more isotropic particles form solids) which has been observed in simulations Strobants\_1 ; Bolhuis ; Veerman  and experiments Maeda  is confirmed by this simple model.

There are two important effects that the restriction of orientations has on the phase diagram topology: (i) The already pointed out stability of the discotic smectic. (ii) The stability of the columnar phase of prolate parallelepipeds for some aspect ratios. The structural properties of inhomogeneous phases that were found through functional minimization allow us to elucidate interesting effects such as the microsegregation behavior of different species in the solids and the depletion effect between particles in the smectics. Those findings endorse the predictive power of the FMF in the description of highly inhomogeneous phases.

###### Acknowledgements.

The author thanks J. A. Cuesta and E. Velasco for useful discussions and a critical reading of the manuscript and E. Jagla for kindly providing his simulations data. This work is part of the research Project No. BFM2004-0180 (DGI) of the Ministerio de Ciencia y Tecnología (Spain). The author was supported by a Ramón y Cajal research contract.


format: text

text:
# On-surface and Subsurface Adsorption of Oxygen on Stepped Ag(210) and Ag(410) Surfaces

## 1 Introduction

The interaction between oxygen and silver surfaces has been the subject of extensive experimental and theoretical research because of its key role in important heterogeneous catalytic reactions, such as partial oxidation of methanol to formaldehyde or ethylene epoxidation . Despite significant efforts, these systems are still not well understood and the identification of the active species involved in these reactions remains unclear. In particular, defect sites and subsurface oxygen species have attracted considerable attention since they are believed to be important in the ethylene epoxidation reaction .

In recent experimental works  Rocca and coworkers have investigated the interaction of oxygen with Ag(410) and Ag(210) surfaces with the aim of understanding the role of steps in the adsorption of oxygen on silver. In particular, using a supersonic molecular beam to dose oxygen on the surface and characterising the final adsorption state by vibrational spectroscopy (HREEL), they have studied the different oxygen species which form on these surfaces. Of particular interest are the results for O/Ag(210) . The HREEL spectra of this system show a peak around 56 meV which is proposed to be due to the vibration of an oxygen atom occupying a subsurface site. Such a peak appears only when preadsorbed oxygen is present on the surface and it is not observed neither on low Miller index surfaces nor on O/Ag(410).

Motivated by these results we have performed an ab initio calculation to investigate the possibility of incorporating oxygen into interstitial subsurface sites on Ag(210) and Ag(410) surfaces.

Recent theoretical works have studied the inclusion of O into subsurface regions of low Miller index silver surfaces, i.e. Ag(111)  and Ag(100) . It was found that pure subsurface O atoms are always unstable with respect to the oxygen molecule but the presence of co-adsorbed on-surface O atoms can stabilise them. More generally, we find that the adsorption in subsurface sites depends significantly on coverage. Isolated subsurface oxygen atoms are always less stable than adatoms occupying favourable on-surface sites, essentially because of the high energy cost for lattice distortion induced by subsurface O incorporation. However, increasing the coverage the preference of on-surface sites decreases and, for example, in O/Ag(111) pure subsurface sites are stabilised at coverages higher than 0.5 ML . This behaviour is due to the fact that the repulsive interaction between the adsorbates is more effectively shielded for the subsurface atoms than for the on-surface ones. In this work we consider high Miller index surfaces in order to investigate the effects of steps on the incorporation of oxygen atoms into subsurface sites.

## 2 Computational Details

Calculations were performed in the framework of density functional theory (DFT) using the generalised gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) . We have used the pseudopotential method with ultra-soft pseudopotentials  and plane-wave basis sets up to a kinetic-energy cutoff of 27 Ry (216 Ry for the charge-density). Details about the Ag and O pseudopotentials are reported in Ref. . Brillouin zone (BZ) integration has been performed with the Gaussian-smearing special-point technique  with a smearing parameter of 0.03 Ry. Calculations have been done using the PWscf package , while the molecular graphics were produced by the XCRYSDEN  graphical package.

We have used periodic super-cells as models of the surface. The Ag(210) surface is modelled with slab of 14 (210) layers, which corresponds to seven (100) layers rotated by an angle $`\varphi _2`$ ($`\mathrm{tan}(\varphi _2)=1/2`$) around the  axis and merged together so as to form the Ag(210) surface. The Ag(410) is modelled with slabs of 20 (410) layers corresponding to five (100) layers rotated by an angle $`\varphi _4`$ ($`\mathrm{tan}(\varphi _4)=1/4`$) and merged together. Adjacent slabs are separated by a vacuum region of at least 16 a.u. Oxygen is adsorbed on both sides of the slab and all the structures have been fully relaxed until the Hellmann-Feynman forces were lower than $`10^3`$ Ry/a.u.

Adsorption of atomic oxygen has been modelled by (2$`\times `$1) surface super-cells. For the Ag(210) surface the BZ integrations have been performed using a (4$`\times `$7$`\times `$1) Monkhorst-Pack mesh, while for the Ag(410) surface a ($`4\times 4\times 1`$) mesh has been applied. All the calculations are non-spin-polarised, because the interaction between O and the Ag substrate results in a negligible spin moment on the O atom.

The chemisorption energies, $`\overline{E}_{\mathrm{chem}}`$, are referred to the clean Ag($`n`$10) surface ($`n=2,4`$) and the isolated oxygen molecule, $`\overline{E}_{\mathrm{chem}}=(E_{\mathrm{O}/\mathrm{Ag}}E_{\mathrm{Ag}}N_\mathrm{O}(E_{\mathrm{O}_2}/2))/N_\mathrm{O}`$, where the total energy of the adsorbate–substrate system, of the clean surface, of the isolated O2 molecule, and the number of adsorbed oxygen atoms are represented by $`E_{\mathrm{O}/\mathrm{Ag}}`$, $`E_{\mathrm{Ag}}`$, $`E_{\mathrm{O}_2}`$ and $`N_\mathrm{O}`$, respectively. $`E_{\mathrm{O}_2}`$ is estimated from a spin-polarised calculation. With this definition, stable adsorbates have negative chemisorption energies.

## 3 Results and Discussion

In previous works on Ag(001) , we found that the surface four-fold hollow site is the most stable one and the chemisorption energy was found to decrease increasing the coverage, a fact related to the repulsive lateral electrostatic interaction between negatively charged oxygen adatoms. Recently, we performed a systematic study of on-surface atomic oxygen adsorption on Ag(410) and Ag(210) surfaces . Figure 1 shows the structure of the two surfaces and the on-surface sites that have been studied. In site A the oxygen adatom lies just above the step and is coordinated with three Ag surface atoms; in site B the oxygen adatom lies in the hollow site just below the step. These two sites are present on both surfaces. The Ag(410) surface possesses also the T1 and T2 sites (see Figure 1) with the O adatom in the four-fold hollow sites on the (100) terrace. We find that the configuration in which the adatoms occupy the sites A forming –O–Ag–O–Ag– rows (configuration labelled A–A) is particularly stable. On the Ag(210) surface the step decoration is significantly more favourable than the other adsorption configurations. On Ag(410), instead, the terrace sites have a chemisorption energy almost degenerate with that of the step sites. In particular, on this surface, when the adatoms are far apart from each other the hollow sites on terrace are slightly more stable than the step edge sites, while at higher coverage the adatoms slightly prefer to decorate the steps.

Let us focus now on the energetics of oxygen inclusion into subsurface interstitial sites, addressing in particular the effects of on-surface co-adsorption. Our results are summarised in Table 1 and discussed in the next paragraphs.

The inclusion of oxygen into subsurface sites is strongly unfavourable with respect to on-surface O adsorption, the resulting chemisorption energies on Ag(210) surface being as low as $`0.07`$ eV for subsurface-octahedral site (Octa) and $`0.14`$ eV for subsurface-tetrahedral site (Tetra).111These calculations were performed with Ag slabs composed of only 10 (210)-layers and somewhat smaller ($`3\times 5\times 1`$) Monkhorst-Pack mesh. It is interesting to observe that the adsorption in Octa and Tetra sites on Ag(100) at $`\theta =0.25`$ ML  is considerably less stable (0.79 eV and 0.69 eV, respectively) than on Ag(210) surface. This indicates that the presence of steps induces a stabilisation of subsurface oxygens. In the right panels of Figure 2 we show top views of the optimised structures with O atom in the Octa and Tetra sites. Note that due to a significant distortion of the substrate lattice induced by subsurface O atoms, the surface structure of Ag(210) can hardly be recognised. However these subsurface sites are significantly stabilised when a sufficient amount of O adatoms is present on the surface. The chemisorption energies for A–A+Octa configuration222 In this configuration O adatoms decorate the step-edge and an additional O atom per ($`2\times 1`$) super-cell is located into Octa site situated straightly below a step-edge Ag atom. are $`0.54`$ and $`0.55`$ eV for Ag(210) and Ag(410) surfaces, respectively. Also observe that the substrate reconstruction is considerably reduced in these cases, and that the Ag(210) and Ag(410) surface structures can be easily recognised from the corresponding optimised O/Ag structures (see left panel of Fig 3 and right panel of Fig 4).

It is interesting to compare the stability of mixed on-surface + subsurface configurations with the purely on-surface ones at a given (high enough) oxygen coverage. In all the configurations we have considered, O adatoms decorate the step-edge (A–A arrangement). An additional oxygen atom is then considered to be located either into subsurface Octa site (A–A+Octa) or adsorbed on a given surface hollow site (see Figs 3 and 4). For Ag(210), we consider only the A–A+B configuration with an additional O adatom located in the on-surface hollow B site situated below the step edge. For Ag(410), in addition to the A–A+B geometry, we also consider the configuration with an additional O atom adsorbed on the terrace T1 site (configuration A–A+T1). On Ag(210) the mixed on-surface + subsurface A–A+Octa configuration is more stable than the pure on-surface A–A+B configuration, $`\overline{E}_{\mathrm{chem}}`$ being $`0.54`$ and $`0.40`$ eV, respectively. On the contrary, the opposite is true for the Ag(410) surface, where the chemisorption energies for the A–A+Octa and A–A+B configurations are $`0.55`$ and $`0.68`$ eV, respectively. This qualitative difference is related to the different O coverages of the two surfaces, which is larger on Ag(210) surface ($`\theta =3/4`$ ML) than on Ag(410) ($`\theta =3/8`$ ML). This suggests that the electrostatic repulsion between the adsorbates in the A–A+B configuration is stronger on Ag(210) than on Ag(410) and makes this geometry on Ag(210) significantly less favourable. Also note that the effect of high coverage in Ag(210) is weaker for the A–A+Octa configuration essentially because the distance between the adatoms is larger than in the A–A+B arrangement and because the interaction involving subsurface O species is shielded more effectively than for purely on-surface O atoms. In addition to this, we have found that on Ag(410) the incorporation of oxygen into subsurface sites is less stable than on-surface adsorption also because of the availability of stable terrace sites. The A–A+T1 configuration is indeed more stable than the A–A+Octa one.

## 4 Conclusions

Our results are compatible with the suggestion, put forward in Ref. , that the HREELS peak found at 56 meV in Ag(210), and missing in Ag(410), is due to atomic oxygen adsorbed subsurface. Indeed our calculations, although performed at different coverages, indicate that subsurface O is not stable in Ag(410), while it is stabilised in Ag(210) at high enough coverage. Further theoretical investigations spanning a larger range of coverages are called for. Lattice vibrational calculations on selected O/Ag($`n`$10) configurations would probably fully clarify the picture.

## Acknowledgements

This work has been supported by INFM (Iniziativa trasversale calcolo parallelo, Sezioni F e G) and by the Italian MIUR (PRIN). All numerical calculations were performed on IBM-SP3 and IBM-SP4 computers at CINECA in Bologna (Italy).


text:
# High Speed Chaos in Optical Feedback System with Flexible Timescales

## I Introduction

Time delay systems are widely used as generators of chaos in applications such as chaos communication  and chaos control . Many of the experiments conducted so far employed lasers with delayed feedback , owing to the simplicity of implementation and feasibility of extension to high-speed operation. Delayed-feedback laser systems also have the potential to generate very high-dimensional and complex chaotic dynamics  and strategies for controlling fast chaos exist for these type of systems .

In this paper we describe a new fast optical time-delay feedback device with flexible dynamical timescales and complexity. Adjusting the pass-band characteristics of the feedback loop allows tuning of the characteristic time scale while the time-delay and the feedback strength control the complexity of the dynamics. The device belongs to the class of optical systems with passive nonlinearity in the feedback loop (see Fig. 1). In our device a Mach-Zehnder interferometer is the source of nonlinearity while the semiconductor laser that provides the optical power acts as a linear current-to-optical-frequency converter. The nonlinearity of the interferometer coupled with the delay in the feedback loop combine to produce a range of steady state, periodic, and chaotic behavior.

The flexible time scale allows us to operate the device at moderate speeds to perform detailed system characterization while it can be also operated at high-speed for applications like chaos communication and control of fast chaos. In this paper we tune the device to a moderate speed so that it generates dynamics with frequencies of hundreds of MHz. Another advantage of the device is that the nonlinearity (interferometer) is easily accessible and tunable so that it can be reproduced and controlled accurately, without involving the internal dynamics of semiconductor lasers as in other devices. Furthermore, this system is constructed with relatively inexpensive components making it an economical choice as a chaotic optical source in future applications.

A distinguishing feature of our device is that it uses an AC-coupled amplifier in the feedback loop so that it has bandpass characteristics. It can thus operate, in principle, in the radio frequency or microwave range using readily available components. This contrasts previous research that used DC-coupled low-pass filter components . Time-delay dynamical systems with bandpass filtering has only recently been investigated and has the advantage that the bandwidth of the chaotic signal can be tailored to fit the desired communication band .

The goal of this paper is to present details about the experimental implementation of our new device and to develop a model that allows us to investigate the nonlinear dynamics of the system. A thorough characterization of the system and a good model are prerequisites for investigating applications such as control of fast chaos, that will be reported elsewhere . We describe the experiment in Sec. II and develop a deterministic model for the device in Sec. III. Subsequently, we investigate in Sec. IV the Andronov-Hopf bifurcation of the steady state. Finally, we will present evidence that our opto-electronic feedback system generates high-dimensional chaos in Sec. V.

## II Experimental setup

In this section we describe details of the experimental implementation of the active laser interferometer with AC-coupled feedback. The device consists of the laser, that acts a current controlled source, the interferometer, that constitutes the passive nonlinearity in the system, and the feedback loop with bandpass characteristics. A schematic of the experimental setup is shown in Fig. 1 where the labels A-L correspond to components that we refer to and describe below.

The light source is an AlGaInP diode laser (A - Hitachi HL6501MG, wavelength 0.65 $`\mu m`$) with a multi-quantum well structure. The diode is housed in a commercial mount (B - Thorlabs TCLDM9) equipped with a bias-T for adding an RF component to the injection current. Thermoelectric coolers in the mount are connected to a proportional-integral-derivative feedback controller (Thorlabs TEC2000) to provide $`1`$ mK temperature stability thereby minimizing frequency and power drift. The output light of the laser is collimated by a lens (D - Thorlabs C230TM-B, f=4.5mm) producing an elliptical beam (1 mm $`\times `$ 5 mm) with a maximum output power of 35 mW.

The passive nonlinearity in the experiment consists of a Mach-Zehnder interferometer with unequal path lengths (path difference 45 cm) into which the laser beam is directed. A silicon photodetector (E - Hammamatsu S4751, DC-750 MHz bandwidth, 15 V reverse bias) measures the intensity of light emitted from one output port of the interferometer. The size of the photodiode is much smaller than the width of the laser beam so only a fraction of the interferometer’s output is detected. The small detector size ensures that only one fringe appears within the beam cross section thus compensating for wavefront aberrations and slight laser beam misalignment and improving the fringe visibility. A neutral density filter is fixed to the front of the laser mount limiting the optical power reaching the photodiode to prevent saturation.

The feedback-loop photodiode produces a current proportional to the optical power falling on its surface. The current flows through a $`50\mathrm{\Omega }`$ resistor. The voltage across that resistor is transmitted down a coaxial cable (F - RU 58, total length $`327`$ cm). The signal emanating from the cable passes through a low-noise, fixed-gain, AC-coupled amplifier (G - MiniCircuits ZFL-1000LN, bandwidth 0.1-1000 MHz), a DC-blocking chip capacitor (H - 220 pF), an AC-coupled amplifier (K - Mini-Circuits ZFL-1000GH, bandwidth 10-1200 MHz), and a second DC-blocking chip capacitor (L - 470 pF). The capacitors reduce the loop gain at frequencies below $`7`$ MHz where a thermal effect enhances the laser’s sensitivity to frequency modulation . The resulting voltage is applied to the bias-T (B) in the laser mount. The bias-T converts the signal into a current and adds it to a DC injection current from a commercial laser driver (C - Thorlabs LDC500).

The entire system is fixed on an optical table using short (2 inch) mounts for mechanical stability. This stability is extremely important as variation in the path length on the order of the wavelength of the laser light (0.65 $`\mu `$m) produces significant power variations at the output of the interferometer. Furthermore, the entire apparatus is covered by an insulating box to prevent thermal expansion or contraction of the mirror mounts due to air currents.

## III Mathematic model of the opto-electronic device

To obtain a model of the device, we consider in turn the relevant physics of the laser diode, the Mach-Zehnder interferometer, and the feedback loop components.

### III-A Source: Semiconductor Laser

The injection current applied to the laser diode is a combination of the DC-bias current and the high-frequency currents due to the time-delayed output of the feedback loop. Modulating semiconductor lasers by varying the input current results primarily in changes of the laser frequency and, to a lesser extend, the laser power. One mechanism relating the input current and frequency shifts is the change of carrier density in the laser device as result of the modulation. A changed carrier density shifts the refractive index of the material that makes up the laser cavity and thereby changes the frequency of the lasing mode. A second mechanism that leads to changes in the optical frequency is a thermal effect that enhances frequency modulations below $`5`$ MHz by as much as a few orders of magnitude

If the pumping current is modulated at a rate significantly slower than the GHz internal time scale typical for semiconductor lasers, then the output power and frequency will adiabatically follow the input so that they depend in a linear fashion on the current when the laser is operated far above threshold.

To determine the main internal time scale of the laser dynamics in the absence of time-delay feedback we determine its relaxation oscillation frequency, $`\mathrm{\Omega }_R`$, by measuring the peak in the power spectrum of the intensity noise using a high-speed spectrum analyzer (Hewlett-Packard 8566B, 22 GHz bandwidth). We find, as expected , that $`\mathrm{\Omega }_R^2ii_{th}`$ where $`i`$ is the DC-injection current and $`i_{th}`$ is the laser threshold injection current. The relaxation oscillation frequency at the nominal operating current of 75 mA used in the experiment is $`\mathrm{\Omega }_R=2.7`$ GHz, and $`i_{th}=40.7`$ mA.

Based on our measurements, we see that the laser will adiabatically follow the injection current for frequencies much less than $`2.5`$ GHz (bounded by $`\mathrm{\Omega }_R`$) and much greater than $`5`$ MHz (bounded by the enhanced frequency tuning due to the thermal effect). Our feedback loop is bandpass limited to within this range and hence the laser can be modeled as a voltage controlled source, a linear device which converts variations of the input voltage $`U(t)`$ that drives the input current into corresponding oscillations of the optical frequency $`\omega (t)`$ and power $`P(t)`$.

The model for the laser we propose is simply

$`P(t)`$ $`=`$ $`\kappa U(t)+P_0,`$ (1)
$`\omega (t)`$ $`=`$ $`\eta U(t)+\omega _0={\displaystyle \frac{\eta }{\kappa }}\left(P(t)P_0\right)+\omega _0,`$ (2)
where $`P_0`$ ($`\omega _0`$) denote the emission power (the optical frequency) of the steady-state, $`U(t)`$ is the voltage applied to the bias-T, and $`\kappa `$ and $`\eta `$ are, respectively, the voltage-to-power and voltage-to-frequency conversion strength.

We note that it is possible to obtain Eqs. (1) and (2) from a direct analysis of the standard rate equations describing semiconductor laser dynamics  in the limit when the amplitude of the injection current modulation is small in comparison to the DC-injection current, which is appropriate for our device. For higher-frequency operation of the device, above $`\mathrm{\Omega }_R`$, the full semiconductor laser rater equations would have to be used in place of Eqs. (1) and (2). Such high-frequency modulation can lead to nonlinear phenomena such as period-doubling cascades, period tripling, and chaos  when the modulation amplitude is large. We do not consider further such behavior.

### III-B Passive Nonlinearity: Mach-Zehnder Interferometer

The passive nonlinearity in our system is an unequal-path Mach-Zehnder interferometer that converts frequency variations of the light emitted by the laser into intensity variations at the photodetector. The photocurrent produced by the detector is converted into a voltage $`V_{det}(t)`$ by a resistance $`r=`$ 50 $`\mathrm{\Omega }`$ and is subsequently amplified. The difference in propagation time between the two paths of the interferometer is $`\mathrm{\Delta }_t=1.5`$ ns, much smaller than the time scale of variation of the laser injection current. Therefore, we can assume that the light waves that reach the detector through the two different paths have the same optical frequency. The voltage $`V_{det}`$ is thus given by

$$V_{det}(t)=\gamma P(t)\left\{1+b\mathrm{cos}\left[\omega (t)\mathrm{\Delta }_t\right]\right\},$$
(3)
where $`\omega (t)`$ is the optical frequency. The parameter $`\gamma `$ describes the overall feedback strength which is determined by the fraction of the power in the laser beam that actually falls onto the detector, the sensitivity of the photodiode, and the electronic amplification of the signal. The parameter $`b[0,1]`$ is the fringe visibility and is defined as

$$b=\frac{P_{max}P_{min}}{P_{max}+P_{min}},$$
(4)
$`b=1`$ corresponds to the ideal situation where the light waves in the interferometer are perfect plane waves and where the beamsplitters divide them exactly in half. In our experiment we find $`b0.8`$.

The optical frequency is related to the observed output power by Eq. (2), which allows us to rewrite Eq. (3) as

$$V_{det}(t)=\gamma P(t)\left\{1+b\mathrm{sin}\left[\alpha \left(P(t)P_0\right)+\varphi \right]\right\},$$
(5)
Here, the parameter $`\alpha =\eta \mathrm{\Delta }_t/\kappa `$ and the constant offset $`\varphi =\omega _0\mathrm{\Delta }_t\pi /2`$ are tunable. By varying the DC component of the injection current, we adjust $`\varphi `$ to $`0`$. The constant parameter $`\alpha `$ determines the sensitivity of the interferometer and can be tuned by varying the path imbalance $`\mathrm{\Delta }_t`$.

The useful range of path-length-differences is limited by the phase noise of the laser light, which is a consequence of the quantum process of spontaneous emission. Increasing the sensitivity of the interferometer by increasing the path difference proportionally increases the effect of phase noise. We chose to set the sensitivity to a value such that the amplitude of noise fluctuations is less than 10% of $`P_0`$, corresponding to $`\alpha `$ 1.9 mW<sup>-1</sup>.

### III-C Electronic Feedback Loop

The electronic feedback loop connecting the detector and diode laser introduces a time-delay. The total roundtrip delay-time $`\tau `$ is somewhat larger than this time-delay, because $`\tau `$ includes, for instance, the free-space propagation time of the laser light. However, for the purpose of modeling, we may assume all processes to be instantaneous and assign the total roundtrip delay-time to the propagation of the signal through the feedback loop. The frequency-limiting effects of the different components in the feedback loop are modeled by a combination of single-pole low- and high-pass filters. The feedback loop is described by

$`\tau _l\dot{V}(t)`$ $`=`$ $`V(t)+V_{det}(t\tau ),`$ (6)
$`\dot{U}(t)`$ $`=`$ $`{\displaystyle \frac{U(t)}{\tau _h}}+\dot{V}(t),`$ (7)
where $`\tau _l`$ ($`\tau _h`$) is the low-pass (high-pass) filter time constant, $`V_{det}`$ is the voltage output of the photodiode, and $`V`$ and $`U`$ are the voltages at the output of the low- and high-pass filters, respectively. The voltage $`U(t)`$ is used to generate the time-dependent current that is injected into the laser diode. Low-pass filtering in the experiment is provided by the limited bandwidth of the photodiode and the electrical connections to the laser. High-pass filtering is due to the two capacitors and the bias-T (see Fig. 1).

Using Eq. (1) we can rewrite Eq. (7) in terms of the laser output power $`P(t)`$

$$\dot{P}(t)=\frac{P(t)P_0}{\tau _h}+\kappa \dot{V}(t).$$
(8)
We determine the parameters $`\tau _l`$, $`\tau _h`$, $`\kappa `$ by measuring the open-loop transfer characteristic of the electronic feedback loop. A signal generator is connected in place of the photodetector so that a sinusoidal voltage $`V_{in}`$ of known amplitude and frequency is injected into the opened feedback loop and eventually into the laser. The optical power is measured directly at the laser output. Figure 2 shows the result of the experiment and the best fit of the model (Eq. (8) and Eq. (6) with $`V_{det}`$ replaced by $`\xi V_{in}`$, where $`\xi `$ is the electronic amplification) with parameter values given in Tab. I. It is seen that the theoretical curve fits to within a few percent of the experimental data everywhere except in the region between 80 and 140 MHz. The dip in the response is caused by either the bias-T or the electronics in the commercial laser mount for which we do not have a detailed circuit diagram.

Aside from this discrepancy, the simple model of the feedback based on single-pole low- and high-pass filters fits the experiment rather well. We choose to ignore the remaining discrepancy because the model successfully reproduces most features of the observed dynamics.

### III-D Full Model and Parameter Values

We obtain a full description of our opto-electronic device by combining Eq. (6) and Eq. (8) with the expression relating the detector voltage to the laser power, Eq. (5). All parameters in this model can be measured and are displayed in Tab.I. The only parameter that was not measured directly is the time-delay $`\tau `$. Direct measurement of this parameter is complicated by contributions of the several electronic components, each of which introduces some unknown phase lag in addition to simple propagation delay. So we use the periodic dynamics of the device to determine $`\tau `$ more precisely. In the next section we present evidence that the steady state becomes unstable through an Andronov-Hopf bifurcation. The frequency of the resulting periodic oscillation depends on the time delay $`\tau `$ and we use the experimentally measured frequency close to the bifurcation point to improve the estimate of $`\tau `$ that we obtain by measuring the propagation delay of the feedback loop.

## IV First Instability: Andronov-Hopf Bifurcation

Our opto-electronic time-delay feedback device can display very complex dynamics. As parameters of the device are changed a series of different bifurcations results in a transition from the initial steady state behavior to chaotic dynamics. In this section we discuss the first such transition in which the steady state becomes unstable and self-sustained periodic oscillations arise.

We know of no exhaustive list of all possible ways that limit cycles (periodic oscillations) can arise in time-delay systems. However, for the well known bifurcations of cycles (those that already exist in two-dimensional systems) this list exists and by examining the scaling of the period and amplitude near the bifurcation one can distinguish between the different bifurcation scenarios . For instance, a supercritical Andronov-Hopf bifurcation is characterized by an amplitude of the stable limit cycle that scales as the square-root of the distance of the bifurcation parameter from the bifurcation point and an oscillation period of finite size that is approximately constant as the bifurcation parameter is varied.

To investigate the bifurcations in our device, we varied the feedback gain $`\gamma `$, which serves as our main bifurcation parameter. A second bifurcation parameter of interest is the delay time $`\tau `$. Experimentally we can change $`\tau `$ by adding or subtracting fixed lengths of coaxial cable to the feedback loop.

For gain values below a critical value $`\gamma <\gamma _C`$, the system is in a steady state with fluctuations of the observed laser output power due only to the inherent phase noise. When the gain is increased through the critical value $`\gamma _C=5.1\pm 0.5`$ mV/mW for $`\tau 19.1`$ ns, the steady state is replaced by periodic oscillations. The dominant frequency of the oscillation is $`51.5\pm 1`$ MHz, which is roughly equal to $`1/\tau `$. This frequency does not change substantially as the gain is increased further. The oscillation amplitude, on the other hand, grows smoothly from zero with increasing gain, as shown in Fig. 3. The spontaneous emission noise of the semiconductor laser leads to an amplification of the amplitude variations (indicated by the larger error bars) close to the bifurcation . It is therefore not possible to pinpoint the bifurcation point exactly and there is no clear $`\sqrt{\gamma \gamma _C}`$ scaling of the amplitude, as expected for an Andronov-Hopf bifurcation. Nevertheless, the smooth amplitude growth and the finite period of the limit cycle at $`\gamma \gamma _C`$ indicate a supercritical Andronov-Hopf bifurcation at $`\gamma _C`$. In the model, which is noise-free, we find an Andronov-Hopf bifurcation at $`\gamma _C=5.34`$ mV/mW.

Next, we experimentally determined the frequency of the limit cycle close to the bifurcation point for different delay times $`\tau `$. In all cases the steady state becomes unstable through an Andronov-Hopf bifurcation. However, we find that the relation $`f1/\tau `$ between the frequency $`f`$ and the delay time $`\tau `$ holds only for a limited range of $`\tau `$. Figure 4 summarizes the relation between $`f`$ and $`\tau `$ that we obtain from experimental (triangles) and numerically calculated (circles) time series. The data suggest that the device transitions from a steady state to limit cycle oscillations with frequencies roughly $`n/\tau `$, where $`n=n(\tau )`$ can be $`1,2,3,\mathrm{}`$.

The origin of the fundamental frequency ($`n=1`$) can be understood intuitively by considering whether a wave circulating in the feedback loop will reinforce itself. If the feedback is positive, a wave will reinforce itself if an integer number of wavelengths equals the propagation length in the loop. If the feedback gain is negative, the propagation length must be a half integer number of wavelengths. In the experiment, we can achieve negative gain by tuning the offset $`\varphi `$ in Eq. (5) such that $`\varphi \pi `$. For this setup we observe, as expected, $`f1/(2\tau )`$ for $`\tau 19.1`$ ns (data not shown).

The appearance of modes with $`n>1`$ for larger delay times $`\tau `$ is due to the fact that the gain in the feedback loop is not perfectly flat over the pass-band. Thus, as the gain is increased from a low level, one particular wavelength will first reach the threshold where the gain in the loop balances the losses. In a system with only low-pass feedback, the gain is highest at low frequencies, so the mode with the lowest frequency is always the one that becomes stable first independent of the delay. On the other hand, the high-pass filter introduces a bias toward high frequencies. For each mode $`n`$ the frequency scales as $`f\tau ^1`$. This implies that the damping effect of the high pass filter on a particular mode becomes more pronounced with increasing delay time $`\tau `$. Finally, a higher order mode $`n>1`$, one that has a higher frequency for a given delay, will reach the stability threshold first. This explains the different modes of periodic oscillations observed in the experiment. In addition, it follows from this argument that there exist specific delay times for which two modes with different frequencies have equal threshold gain $`\gamma _C`$ (double Hopf point). However, since none of the delay-times used in the experiment are close to one of these special $`\tau `$ we do not consider further this case.

## V Chaos

Beyond the Hopf bifurcation successively more complex dynamics develops as the gain is increased, as shown in Figure 5. At feedback gains higher than the Andronov-Hopf bifurcation point, the initially sinusoidal oscillations begin to square off, as shown in Fig. 5a. The square shape of the waveform results in prominent odd harmonics in the spectrum (Fig. 5b). As the gain increases, a small, broad peak appears at about half the fundamental frequency as shown in Fig. 5d. The peak at roughly half the fundamental frequency is three orders of magnitude below the fundamental. The weakness and broadness of this peak coupled with the presence of phase noise may explain why no clearly period doubled behavior is apparent in the time domain (Fig. 5c). As the gain is increased further, the broad background rises and the tall peaks at the fundamental frequency and its harmonics weaken. The power spectrum for $`\gamma =17.6`$ mV/mW, shown in Fig. 5f, is quite broad and the peaks have nearly dropped to the level of background which has risen significantly above the noise floor ($`2\times 10^3`$ mW<sup>2</sup>/MHz<sup>2</sup> measured just below the Andronov-Hopf bifurcation point). This is indicative of high dimensional chaos in the system.

A similar very broad and featureless power spectrum in the chaotic regime for an optical system with passive nonlinearity and bandpass feedback was reported in Ref. . There, the authors synchronize two devices and successfully communicate information, thus demonstrating that the cause of the broadband spectrum is deterministic chaos. Because of the similarity of their device to ours, we believe that the observed complex behavior for large gain values shown in Fig. 5e is due to chaotic deterministic dynamics.

To support this claim we show in Fig. 6 time series and power spectra obtained by numerical simulation of the deterministic model<sup>1</sup><sup>1</sup>1Time series from the model are obtained using an Adams-Bashforth-Moulton predictor-corrector algorithm.. The match with the experimental data is good, as can be seen by comparing Fig. 5 to Fig. 6, despite the simplicity of the model and the uncertainty in the estimation of the parameters. Figure 6 also shows Poincaré sections obtained by recording the location where the trajectory uni-directionally crosses the plane $`V(t)=\gamma P_0`$ in the three-dimensional space spanned by $`(V(t),P(t),P(t\mathrm{\Delta }))`$ with $`\mathrm{\Delta }<\tau `$. The simulations confirm that the system is on a limit cycle for $`\gamma =9.4`$ mV/mW, which is clear from the power spectrum (Fig. 6b) and immediately obvious in the Poincaré section (Fig. 6c). They also show that the limit cycle bifurcates to a torus-attractor for increased gain ($`\gamma =13.2`$ mV/mW) appearing as closed curve in the Poincaré section (Fig. 6f). The power spectrum (Fig. 6e) exhibits a comb-like structure due to the two incommensurate frequencies of the quasi-periodic oscillation. Note, that there is not only a strong peak at $``$26.6 MHz, roughly half the fundamental frequency, but a definite peak at 1.8 MHz. This is well below the 3 dB cutoff point of the high-pass filter. At present we do not understand the origin of this low-frequency feature and cannot explain why the main new frequency component should be so close to one half the fundamental frequency. For even larger gains the system is chaotic, with a very broadband spectrum (Fig. 6h) and no discernible structure in the Poincaré section (Fig. 6i). We computed the spectrum of Lyapunov exponents for the model with $`\gamma =17.6`$ mV/mW and find that the largest few exponents are positive (details on the computation of the Lyapunov spectrum can be found in Ref. ). This proves that the model dynamics shown in Fig. 6g are on a chaotic attractor. We estimate the attractor’s Lyapunov dimension to be $`D_L22`$.

In this section we have presented evidence that our device exhibits high dimensional chaos. We showed that the device undergoes a sequence of bifurcations from steady state to aperiodic oscillations with a broad and featureless power spectrum. We were able to reproduce this behavior with the deterministic model and we quantified the dimensionality of the chaotic attractor.

## VI Summary and Discussion

We describe an optical feedback device that can produce high-dimensional chaos and that allows adjustment of the characteristic time scales of the oscillations by changing the bandpass characteristics in the feedback loop. The nonlinearity in the device is accessible and reproducible. We develop a simple model that allows quantitative predictions about the behavior of the physical device and use it to determine the critical gain and frequency of the Andronov-Hopf bifurcation of the steady state. We observe that the device transitions to chaos with a very broadband frequency spectrum and find that this matches the model behavior.

We find that the inclusion of a high-pass filter significantly changes the qualitative dynamics of optical feedback systems with passive nonlinearity in comparison to only low-pass filtering as in the Ikeda system . Bandpass feedback allows not only “fundamental” frequencies $`f(2\tau )^1`$ but oscillations with $`f\tau ^1`$ become possible. The route to chaos is apparently changed when the feedback of DC-signals is blocked. That is, we do not observe a period doubling route to chaos but a more complicated transition, the details of which are not yet fully understood.

This chaotic opto-electronic device is ideally suited for both experimental investigation of fast nonlinear dynamics and technological application of high-speed chaos. For example, we use it to investigate control of fast chaos and are able to successfully stabilize a periodic orbit with a period of 12 ns, faster than any reported . This work will be described in a later paper . Also, by adjusting the time-scale of the oscillation the device could be made to oscillate at GHz frequencies suitable for use in a practical chaos communication system.


text:
# A scheme for dense coding in the non-symmetric quantum channel

## I Introduction

The quantum entanglement state  is an important tool distinguishing the quantum mechanics from the classical physics. Due to the non-local properties of quantum entanglement state, it has been widely used in many fields such as quantum teleportation s3 , quantum dense coding s4  and quantum key distribution s5 . Here we will focus our attention on quantum dense coding. In 1992, Bennett and Wiesner proposed the first scheme for the quantum dense coding s4 , in which Bell states are used as quantum channel. In 1996, quantum dense coding was experimentally presented by Mattle $`et`$ $`al`$ in an optical system s6 . Recently Liu $`et`$ $`al`$ presented a protocol for dense coding with multi-level entangled states s7 . However in these schemes, the quantum channels are symmetric, that is to say, the dimension of the Hilbert space of the particle with sender is the same as that of the particle with receiver. In 1998, Bose $`et`$ $`al`$ generalized Bennett and Wiesner’s scheme for dense coding into multi-parties, one receiver intends to receive messages from $`N`$ senders s8  . In this case, the quantum channel is non-symmetric. The particles with senders are in $`2^N`$\- dimensional Hilbert space, while the particle with 2-dimensional Hilbert space belongs to receiver. In present paper, we give a scheme for quantum dense coding in the case of non-symmetric Hilbert spaces of the sender and receiver’s particles sharing the quantum maximally entangled state, and consider the efficiency of classical information gain.

## II A scheme for dense coding in the non-symmetric quantum channel

In Bennett and Wiesner’s first scheme for the quantum dense coding, sender Alice and receiver Bob share a pair of entangled particles in the Bell state. Alice performs one of the four 1-qubit unitary operations given by the identity $`I`$ or the Pauli matrices $`(\sigma _x,i\sigma _y,\sigma _z)`$ on her particle. Each of the unitary operations maps the Bell state to a different member of the four Bell states. Then Alice sends her particle to Bob. Bob can obtain two bits of classical information from the joint measurement on two particles. Now we would like to generate the Bennett and Wiesner’s symmetric scheme to the non-symmetric one.

To present our scheme clearly, let us first begin with dense coding between two parties in $`3\times 2`$-dimension. Suppose the particle 1 in 3-dimensional Hilbert space belongs to Alice, and the particle 2 in 2-dimensional Hilbert space is with Bob. They share the maximally entangled state:

$$|\mathrm{\Psi }_{00}_{12}=\frac{1}{\sqrt{2}}(|00+|11)_{12}.$$
(1)
Through simple calculation, it can be shown that the single-body operators on particle 1:

$$\begin{array}{c}U_{00}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],U_{01}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],\hfill \\ U_{10}=\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right],U_{11}=\left[\begin{array}{ccc}0& 0& 1\\ 1& 0& 0\\ 0& 1& 0\end{array}\right],\hfill \\ U_{20}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right],U_{21}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right]\hfill \end{array}$$
(2)
will transform $`|\mathrm{\Psi }_{00}`$ into the corresponding state respectively:

$$U_{00}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|00+|11)=|\mathrm{\Psi }_{00},$$
(3)
$$U_{01}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|00|11)=|\mathrm{\Psi }_{01},$$
(4)
$$U_{10}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|10+|21)=|\mathrm{\Psi }_{10},$$
(5)
$$U_{11}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|10|21)=|\mathrm{\Psi }_{11},$$
(6)
$$U_{20}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|20+|01)=|\mathrm{\Psi }_{20},$$
(7)
$$U_{21}|\mathrm{\Psi }_{00}=\frac{1}{\sqrt{2}}(|20|01)=|\mathrm{\Psi }_{21}.$$
(8)
Alice operates one of the above unitary transformations, and sends her particle to Bob. Bob takes only one measurement in the base $`\{|\mathrm{\Psi }_{00},|\mathrm{\Psi }_{01},\mathrm{},|\mathrm{\Psi }_{21}\}`$, and he will know what operation Alice has done, or to say , what the messages Alice has encoded in the quantum state. As a result, Bob gets $`\mathrm{log}_26`$ bits of information through only one measurement. Thus the dense coding is realized.

It is straightforward to generalize the above protocol to arbitrarily different dimension for two parties. Assume that the Hilbert space of Alice’s particle 1 is non-symmetric with that of Bob’s particle 2, i.e. the dimension $`p`$ of the Hilbert space of particle 1 is not the same as the dimension $`q`$ of particle 2. Without loss of generality, we choose $`p>q`$. The non-symmetric quantum channel shared by Alice and Bob is the maximally entangled state:

$$|\mathrm{\Psi }_{00}_{12}=\frac{1}{\sqrt{q}}(|00+|11+\mathrm{}+|q1q1)_{12}.$$
(9)
Clearly, an orthogonal base can be composed of the general Bell states of particle 1 and 2

$$|\mathrm{\Psi }_{mn}=\underset{j}{}\mathrm{e}^{2\pi \mathrm{i}jn/q}|(jm)\mathrm{mod}p_1|j_2/\sqrt{q},$$
(10)
where $`m=0,1,\mathrm{},p1;`$ $`n,j=0,1,\mathrm{},q1`$. The single-body operators on particle 1 satisfying $`U_{mn}|\mathrm{\Psi }_{00}=|\mathrm{\Psi }_{mn}`$ may be explicitly written out:

$$U_{mn}=\underset{j}{}\mathrm{e}^{2\pi \mathrm{i}jn/q}|(jm)\mathrm{mod}pj|,$$
(11)
where $`m=0,1,\mathrm{},p1;`$ $`n,j=0,1,\mathrm{},q1`$. The procedure for realizing dense coding in this case is similar to the above one. In order to send information, Alice performs one of the operations in (11) to her particle and sends it to Bob. Bob then performs a collective measurement in the general Bell base $`\{|\mathrm{\Psi }_{mn}_{12}|m=0,1,\mathrm{},p1;q=0,1,\mathrm{},q1\}`$ to find out what Alice has done to particle 1 and therefore reads out the encoded message. In this case, Bob can gain $`\mathrm{log}_2p\times q`$ bits of information. If $`p=q`$, the quantum channel becomes symmetric, correspondingly the bits of information become $`\mathrm{log}_2p^2=2\mathrm{log}_2p`$.

## III improving the efficiency of the dense coding by entanglement transfer

When the quantum channel between Alice and Bob is symmetric, besides a more level particle is with Alice, then the quantum channel may be transformed into non-symmetric one, and hence improve the efficiency of the dense coding. For the sake of clearness, we start our discussion from the following simple case.

Suppose the initial quantum channel between Alice and Bob is

$$|\mathrm{\Phi }_{12}=\frac{1}{\sqrt{2}}(|00+|11)_{12},$$
(12)
the particle 1 belongs to Alice, and the particle 2 to Bob, they are both in 2-dimensional Hilbert space. Besides the particle 3 belonging to Alice is in 3-dimensional Hilbert space, its initial state is $`|0_3`$. So the initial state of the three particles is

$$|\mathrm{\Psi }_{123}=\frac{1}{\sqrt{2}}(|00+|11)_{12}|0_3.$$
(13)
First Alice applies the following unitary operation on her particles 1 and 3 under the base $`\{|00_{13},|01_{13},|02_{13},|10_{13},|11_{13},|12_{13}\}`$:

$$U_{13}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$$
(14)
correspondingly $`|\mathrm{\Psi }_{123}`$ becomes

$$U_{13}|\mathrm{\Psi }_{123}=\frac{1}{\sqrt{2}}(|00_{13}|0_2+|11_{13}|1_2).$$
(15)
Second Alice applies the following unitary operation on her particles 1 and 3:

$$U_{13}^{}=\left[\begin{array}{cccccc}1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right],$$
(16)
and hence equation (15) becomes

$$\begin{array}{c}U_{13}^{}\frac{1}{\sqrt{2}}(|00_{13}|0_2+|11_{13}|1_2)\hfill \\ =\frac{1}{\sqrt{2}}(|00_{23}+|11_{23})|0_1\hfill \\ =|0_1|\mathrm{\Phi }^{}_{23}.\hfill \end{array}$$
(17)
As a result, the entanglement between the particles 1 and 2 is transferred to the entanglement between the particles 2 and 3. Before this, using the quantum channel between the particles 1 and 2 for dense coding, Bob gets $`\mathrm{log}_24`$ bits of information, but after that, Bob can gain $`\mathrm{log}_26`$ bits of information, so the efficiency of information transmission is increased.

In the following section, we consider the general case. Assume that Alice and Bob share the maximally entangled state:

$$|\mathrm{\Phi }_{00}_{12}=\frac{1}{\sqrt{q}}(|00+|11+\mathrm{}+|q1q1)_{12},$$
(18)
where the particle 1 is with Alice, and the particle 2 with Bob, they are both in $`q`$-dimensional Hilbert space. Besides the particle 3 belonging to Alice is in $`p`$-dimensional Hilbert space, its initial state is $`|0_3`$, and $`p>q`$. Therefore three particles will in the state $`|\mathrm{\Phi }_{00}_{12}|0_3`$. Alice performs the following unitary operations on her particles 1 and 3:

$$U_{13}=\underset{i=0}{\overset{q1}{}}|iii0|+\underset{j=1}{\overset{q1}{}}|j0jj|+\underset{m,n}{}|mnmn|,$$
(19)
where $`m=0,1,\mathrm{},q1`$; $`n=1,2,\mathrm{},p1,`$ and $`mn`$;

$$U_{13}^{}=\underset{i=0}{\overset{q1}{}}|0iii|+\underset{j=1}{\overset{q1}{}}|jj0j|+\underset{m,n}{}|mnmn|,$$
(20)
in the sum of above equation, when $`m=0,`$ we let $`n=q,q+1,\mathrm{},p1`$; when $`m=1,2,\mathrm{},q1`$, then $`n`$ should take $`n=0,1,\mathrm{},p1`$ with the condition $`mn`$. After that $`|\mathrm{\Phi }_{00}_{12}|0_3`$ correspondingly becomes

$$\begin{array}{c}U_{13}^{}U_{13}|\mathrm{\Phi }_{00}_{12}|0_3\hfill \\ =\frac{1}{\sqrt{q}}(|00+|11+\mathrm{}+|q1q1)_{23}|0_1.\hfill \end{array}$$
(21)
So the entanglement transfer between the particles 1 and 2 and the particles 2 and 3 has been successfully realized. In other words, a symmetric quantum channel has been transformed into the non-symmetric one. Clearly, before this, Bob gets $`\mathrm{log}_2q^2`$ bits of information, now Bob can gain $`\mathrm{log}_2pq`$ bits of information, and hence the efficiency of information transmission can be increased.

## IV conclusion and discussion

In summary, we have given a scheme to transform the symmetric level quantum maximally entangled state between two parties into non-symmetric one, and using it as the quantum channel for dense coding. We have also considered the efficiency of information transmission. Comparing our scheme with one for symmetric dense coding, obviously ours is more efficient. We hope this scheme can be realized by the experiment.

###### Acknowledgements.

This work was supported by Hebei Natural Science Foundation under Grant No. 101094.