text:
# REVTeX 4 Prediction of Neutron Elastic Form Factors Using GPDs from Proton Elastic Form Factors and Isospin Symmetry

(June 29, 2003)

## Abstract

The elastic neutron form factors $`G_{En}`$ and $`G_{Mn}`$ are calculated in a GPD framework using GPDs obtained from fits to proton elastic form factors $`G_{Ep}`$ and $`G_{Mp}`$, and isospin symmetry, with no further changes in parameters. The results for $`G_{En}`$ are in good agreement with existing data, while those for $`G_{Mn}`$ are fair. The calculations predict the form factors for future measurements at higher $`Q^2`$.

form factors; generalized parton distributions; baryon structure

In recent years, the development of generalized parton distributions (GPDs) ji ; rad\_gpd ; collins  has opened the possibility of describing a great variety of exclusive reaction in the multi GeV range in terms of a common nucleon structure. The constraints imposed by the description of many types of reactions offers the possibility of modeling the longitudinal and transverse parton structure of nucleons.

Among the most direct consequences of the GPD formalism are the sum rules which relate the various GPDs to the hadronic form factors. Thus the proton elastic helicity conserving and helicity-flip form factors may be written, respectively, as :

$$F_{1p}(t)=_1^1\underset{q}{}e_qH_p^q(x,\xi ,t)dx$$

$$F_{2p}(t)=_1^1\underset{q}{}e_qE_p^q(x,\xi ,t)dx$$

where $`t=Q^2`$ is the momentum transfer to the proton, $`\xi `$ is the longitudinal momentum transfer, and $`q`$ signifies quark flavors. Without loss of generality one may work in a coordinate system in which the momentum transfer $`t`$ is transverse so that $`\xi =0`$, and the GPDs may be written:

$$H_p^q(x,t)H_p^q(x,t,\xi =0)E_p^q(x,t)E_p^q(x,t,\xi =0)$$

Several authors rad\_wacs ; kroll ; kroll\_GPD ; burkardt  have modeled the GPDs by Gaussian functions which embody general expected properties. In particular, $`H_p^q(E,t=0)f_p^q(x)`$, the unpolarized quark distribution function and asymptotically $`H^q(E,t\mathrm{})`$ narrows toward $`x=1`$ (see burkardt ; stoler\_gpd ). In terms of a Gaussian a simple model is,

$$H_p^q(x,t)=f_p^q(x)e^{\overline{x}t/4x\lambda _H^2}$$
(1)
in which $`\overline{x}1x`$. For $`E_p^q(x,t)`$ the we take the simple ansatz. afanasev

$$E_p^q(x,t)=k_p^q(x)e^{\overline{x}t/4x\lambda _E^2}.$$
(2)
To account for hard components of $`F_{1p}`$ at $`t>10`$ ref. stoler\_gpd  modified the specific functional form for $`H_p^q(x,t)`$ and $`E_p^q(x,t)`$ as a Gaussian plus small power law shape in $`t`$. 111As in ref. stoler\_gpd , eq.17, the parameter for the power law part of the GPD is $`A_h=0.18`$. Essentially no power law shape was required to fit $`F_{2p}`$, where the data only goes up to $`t6`$ GeV2/c2.

$`H_p^q(x,t)=f_p^q(x)exp(\overline{x}t/4x\lambda _H^2)+\mathrm{}`$ (3)
$`E_p^q(x,t)=k_p^q(x)exp(\overline{x}t/4x\lambda _E^2)+\mathrm{},`$ (4)
in which $`\mathrm{}`$ indicates the addition of small power law components in $`t`$.

To obtain $`E_p^u`$ and $`E_p^d`$, needed for eqs. 1 to 4 , the available data for $`G_{Mp}`$ and the recent JLab data jones ; gayou  on $`G_{Ep}/G_{Mp}`$ were fit, as in ref. stoler\_gpd .

The conditions at $`t`$=0 were also required, i.e. $`H_p(x,0)=e_uf_u(x)+e_df_d(x)`$ and $`E_p(x,0)=k_p^u(x)+k_p^d(x)`$. The valence quark distribution functions $`f_p^u(x)`$ and $`f_p^d(x)`$ are measured in DIS, and obtained from refs. rad\_wacs ; martin . The functions $`k_p^u(x)`$ and $`k_p^d(x)`$ are not obtainable from evaluations of DIS . Following ref. stoler\_delta  the simple phenomenological assumption $`k_p^q(x)\sqrt{1x}f_p^q(x)`$ was used. This results in a satisfactory ratio of $`F_{2p}/F_{1p}`$, since for large $`t`$, the quantity $`\sqrt{1x}1/\sqrt{t}=1/Q`$ with normalization obtained by requiring the proton $`F_{2p}(0)=1.79`$.

Adequate fits to the measured $`G_{Mp}`$ and $`G_{Ep}/G_{Mp}`$, or equivalently $`F_{1p}`$ and $`F_{2p}/F_{1p}`$, were obtained with $`\lambda _H=0.76`$ GeV/c and $`\lambda _E=0.67`$ GeV/c. The results are shown in figs. 1 and 2.

This gives

$$F_{1p}^u(0)=e_uH_p^u(x,0)𝑑x=e_uf_p^u(x)𝑑x=4/3$$

$$F_{1p}^d(0)=e_uH_p^d(x,0)𝑑x=e_df_p^d(x)𝑑x=1/3$$

and

$$F_{2p}^u(0)\kappa _p^u=e_uE_p^u(x,0)𝑑x=k_p^u(x)𝑑x=1.67$$

$$F_{2p}^d(0)\kappa _p^d=e_uE_p^d(x,0)𝑑x=k_p^d(x)𝑑x=2.03$$

The neutron form factors were obtained from the fit to proton form factors by applying isospin symmetry.

$$H_n^u(x,t)=H_p^d(x,t)$$

$$H_n^d(x,t)=H_p^u(x,t)$$

$$E_n^u(x,t)=E_p^d(x,t)$$

$$E_n^d(x,t)=E_p^u(x,t)$$

and

$$F_{1n}(t)=_1^1\underset{q}{}e_qH_n^q(x,t)dx$$

$$F_{2n}(t)=_1^1\underset{q}{}e_qE_n^q(x,t)dx$$

$$G_{En}=F_{1n}\kappa _n\tau F_{2n}$$

$$G_{Mn}=F_{1n}+\kappa _nF_{2n}$$

with $`\kappa _n`$ = -1.91 $`\mu _N`$.

The result for $`G_{En}`$ is shown in Fig. 3. The calculated form factor is somewhat lower than the existing data in the region $`Q^2=t<0.75`$ GeV2/c2, but accounts well for the new JLab Hall C data for $`Q^2>0.75`$ zhu ; madey . There is excellent agreement with the results of a calculation of ref. miller , which is also shown in the figure. The calculation of ref. miller  uses a completely different framework, employing a relativistic constituent quark model with a pion cloud. For $`G_{En}`$ the pion cloud is important at small $`Q^2`$, where the constituent quark contribution is very small. However, for $`Q^2>1`$ GeV2/c2 the quarks become most important, with the role of the pion cloud diminishing. In the present calculation, the contribution of the sea quark pairs, which presumably would mimic the pion cloud, was set to zero. The importance of a rigorously relativistic calculation of both the constituent quarks and pion cloud is stressed in ref. miller . For example, at high $`Q^2`$ the lower components of the Dirac spinors, which introduce orbital angular momentum, become important. The calculation of ref. miller  employs several parameters, however the $`Q^2`$ dependence of the form factor at higher $`Q^2`$ appears to be governed more by relativistic effects than the specific parameter set used. In particular a large number of sets of these parameters can be found to give similar $`Q^2`$ dependence.

As seen in Fig. 3 both calculations give results at high $`Q^2`$ which lie above the Galster parameterization galster , as do the most recent experimental data madey . This is not surprising since the Galster parameterization is simply an ad hoc fit to low $`Q^2`$ data.

The result for $`G_{Mn}`$ are shown in Fig. 4. Here, the fit to the experimental data is somewhat poorer than for $`G_{En}`$. Also shown is the result of the calculation of ref. miller . Curves are shown for two of the many parameter sets which fit the data. A possible reason for the better fit may be that ref. miller  chooses parameters in such a way that requires the fit to be rather good for all four elastic form factors, while in the present case only the proton form factors are fit, and then isospin symmetry is applied to obtain the neutron form factors with the parameters fixed.

This note has pointed out the usefulness of GPDs in describing elastic form factors. Alternatively, the elastic form factors, together with isospin symmetry can be very important for constraining nucleon structure through GPDs. Further constraints of details of nucleon structure will be possible by including other high $`t`$ experiments into the fit procedure. These include high $`W`$ high $`t`$ real and virtual Compton scattering, and single meson photo and electroproduction, such as described in refs. rad\_wacs ; kroll ; kroll\_huang . It would be quite interesting if conceptual connections could be made between this technique and those of recent relativistic constituent quark models with a pion cloud such as in ref. miller , or recent helicity non-conserving pQCD based approaches ref. ralston ; belitsky  which have had some success in explaining the $`Q^2`$ dependence of $`F_{2p}`$.

Acknowledgments: The work was partially supported by the National Science Foundation.


format: text

text:
# Atomic scale friction between clean graphite surfaces

## Abstract

We investigate atomic scale friction between clean graphite surfaces by using molecular dynamics. The simulation reproduces atomic scale stick-slip motion and low frictional coefficient, both of which are observed in experiments using frictional force microscope. It is made clear that the microscopic origin of low frictional coefficients of graphite lies on the honeycomb structure in each layer, not only on the weak interlayer interaction as believed so far.

Atomic scale friction, Graphite, Molecular dynamics simulation, Frictional force microscope  preprint: \***

Friction is one of the most familiar physical phenomena and has been investigated from long ago due to its importance in various kinds of machinery and for the understanding of dynamics of many systems in science. Last two decades, atomic scale friction has been attracted much attention for the understanding of fundamental mechanisms of macroscopic friction and for many fields of high precision engineering such as nanomachine(Persson, ). Frictional force microscope (FFM) has played an important role in study of atomic scale frictionPersson ; Mate . Behavior of frictional force in FFM experiments has been studied theoretically and numerically based on the model which consists of a single atom tip and the potential by a substrate surfaceFujisawa1 ; Morita ; Sasaki1 ; Sasaki2 ; Hols . The magnitude of load in FFM experiments in layered materials, however, is much higher than that estimated numerically and theoreticallyMate ; Sasaki3 . Some groups pointed out that a flake cleaved from the substrate plays a crucial role in FFM experiments in such systemsMate ; Morita ; Miura1 ; Miura2 ; Miura3 . It is not appropriate to discuss frictional phenomena with a flake by using the single atom tip model. We investigate in this letter atomic scale friction between clean graphite surfaces by molecular dynamics (MD). The model employed here simulates friction in the FFM experiments on a highly oriented pyrolytic graphite substrate with a flake. The present work enables us to understand results of the FFM experiments more deeply. Graphite is one of the most important solid lubricants. It turns out that the microscopic origin of lubrication property of graphite lies on the honeycomb structure in each layer, not only on the weak interlayer interaction.

The model consists of a monolayer graphite substrate, a monolayer graphite flake and a spring which drives the flake as shown in fig. 1(a). Atoms in the flake obey the following eqs. of motion,

$`m_\mathrm{C}\ddot{x}_i^\alpha =\gamma \dot{x}_i^\alpha {\displaystyle \frac{V_{\mathrm{Sub}}}{x_i^\alpha }}{\displaystyle \frac{V_\mathrm{I}}{x_i^\alpha }}{\displaystyle \frac{k_\alpha }{N_\mathrm{f}}}\left(x_\mathrm{C}^\alpha x_\mathrm{B}^\alpha \right).`$ (1)
Here $`m_\mathrm{C}`$ is the mass of the carbon atom, $`x_i^\alpha `$, $`\alpha `$ $`=`$ $`x`$, $`y`$, $`z`$, are components of the $`i`$-th atom position vector and a dot above $`x_i^\alpha `$ denotes the time derivative. The $`x`$ and $`y`$-directions are shown in fig. 1(b) and the $`z`$-direction is normal to the surface plane of the substrate. Atoms in the substrate are fixed because their deformations are negligibly small in the present parameter range. The first term in the r. h. s. of eq. (1) is a damping term proportional to the velocity, which reproduces the energy dissipation to the substrate. We set a damping constant $`\gamma `$ $`=`$ 1.14 $`\times `$ 10-4 nNs/m. This value of $`\gamma `$ holds the stability of the flake for the time step of MD simulation $`\mathrm{\Delta }t`$ = 2.74 $`\times 10^{15}`$ sec. The second term denotes the force coming from interatomic interactions between the flake and the substrate. $`V_{\mathrm{Sub}}`$ denotes the sum of pair atomic interaction potential between an atom in the flake and that in the substrate. We employ the Lennard-Jones $`(12,6)`$ potential with a cut off length $`R`$, outside of which the interatomic potential is negligible. We set Lennard-Jones parameters and the cut-off length as $`\epsilon `$ $`=`$ 0.965 $`\times `$ 10-2 eV, $`\sigma `$ $`=`$ 0.340 nm and $`R`$ $`=`$ 2.0 nm. The above values of $`\epsilon `$ and $`\sigma `$ are often employed in the interlayer atomic potential in bulk graphite Taber . The third term denotes the force by the intraflake atomic interactions. We employ linear spring potentials for the bond length and the angle between bonds of nearest neighbor carbon atoms and $`z`$-displacement of an atom against nearest atoms. They are employed in the study of lattice vibrations and specific heat of graphite (Yoshimori, ) and FFM images (Sasaki1, ; Sasaki2, ). The fourth term denotes the force by the driving spring with the components of spring constants, $`k_\alpha `$, $`\alpha `$ = $`x`$, $`y`$, $`z`$. The spring corresponds to a tip and a cantilever in FFM. $`N_\mathrm{f}`$, $`x_\mathrm{C}^\alpha `$ and $`x_\mathrm{B}^\alpha `$ denote the number of carbon atoms in the flake, $`\alpha `$-components of a center position of the flake mass and a base position of the driving spring, respectively. Here we set $`k_x`$ = $`k_y`$ = 1.5 eV/nm2, $`k_z`$ = 5.0 eV/nm2. Flake size is not known experimentally. The sizes of the flake $`N_\mathrm{f}`$ in the present calculation are 6, 20, 42 and 110. The results shown below are those for $`N_\mathrm{f}=110`$size . The shape of the flake is shown in fig. 1(b). The size dependence is discussed later. $`\stackrel{}{𝒙}_\mathrm{B}`$ obeys the following eq. of motion,

$`R_\mathrm{C}m_\mathrm{C}\ddot{x}_\mathrm{B}^z=\gamma \dot{x}_\mathrm{B}^zL+k_z\left(x_\mathrm{C}^zx_\mathrm{B}^z\right).`$
$`\dot{x}_\mathrm{B}^x=V_\mathrm{S}^x,\dot{x}_\mathrm{B}^y=V_\mathrm{S}^y,`$ (2)
Here $`R_\mathrm{C}`$, the ratio of the effective mass of the tip and the cantilever to the mass of the carbon atom, is assumed to be 100.0. $`L`$ and $`(V_\mathrm{S}^x,V_\mathrm{S}^y)`$ denote the loading force and the lateral driving velocity, respectively. The values of these parameters are kept constant during driving.

Equations of motion, (1) and (2) are solved numerically by using the Runge-Kutta method. In the initial state the system is stable with the applied loading force $`L`$ and no lateral expansions of the driving spring. The configuration corresponds to the AB stacking structure of bulk graphite as shown in fig. 1(b). The frictional force is defined as a driving direction component of the force of the driving spring, $`k_\alpha \left(x_\mathrm{B}^\alpha x_\mathrm{C}^\alpha \right)`$.

It is to be noted that deformations of the flake are negligibly small in the present calculation. This is because that the interlayer atomic interactions between the flake and the substrate under load are much weaker than the intralayer atomic interactions in the flake. It is also to be noted that the rotation of the flake around $`\stackrel{}{𝒙}_\mathrm{C}`$ in the $`x`$-$`y`$ plane is inhibited due to the potential barrier against it. Due to these two features the dynamics of the center of the flake mass, $`\stackrel{}{𝒙}_\mathrm{C}`$, governs the dynamics of the system. At first we focus on the typical motions of $`\stackrel{}{𝒙}_\mathrm{C}`$. Figure 2 shows the driving direction component of the center of the flake mass, $`x_\mathrm{C}`$, as a function of that of the base position of the driving spring, $`x_\mathrm{B}`$, for $`L`$ $`=`$ 100 nN, $`V_\mathrm{S}`$ $`=`$ 0.22 m/s in the case driven along the $`x`$ (a) and $`y`$ (b) directions, respectively. The stick-slip motions of $`\stackrel{}{𝒙}_\mathrm{C}`$ are observed. These periods agree with the periods of the graphite lattice in each direction, 0.426 and 0.246 nm, respectively.

Figure 3 shows trajectories of $`\stackrel{}{𝒙}_\mathrm{C}`$ in the $`x`$-$`y`$ plane for $`L`$ $`=`$ 100 nN, $`V_\mathrm{S}`$ $`=`$ 0.22 m/s. The contour lines denote equipotential lines of the substrate potential. The motion of $`\stackrel{}{𝒙}_\mathrm{C}`$ starts from the initial position A. The solid line in the figure denotes the trajectory of $`\stackrel{}{𝒙}_\mathrm{C}`$ driven along the $`x`$-direction.

By applying driving force $`\stackrel{}{𝒙}_\mathrm{C}`$ goes slowly to the position B, which is the point of inflection of the total potential. As soon as $`\stackrel{}{𝒙}_\mathrm{C}`$ arrives at B, $`\stackrel{}{𝒙}_\mathrm{C}`$ slips to the next sticking position close to D through the position C, both of which are the positions of the substrate potential minimum. The sticking time around D is much shorter than that around A due to large elongation of the driving spring along the $`y`$-direction and much longer than the slipping time. The sticking times around A and D are about 99.88 and 0.08 % of the total period of the stick-slip motion, respectively. As soon as $`\stackrel{}{𝒙}_\mathrm{C}`$ arrives at E, which is also the point of inflection, $`\stackrel{}{𝒙}_\mathrm{C}`$ slips again close to I, which is equivalent to A. Then $`\stackrel{}{𝒙}_\mathrm{C}`$ repeats the periodic stick-slip motion. The dashed line shown in the figure denotes the trajectory of $`\stackrel{}{𝒙}_\mathrm{C}`$ in the case driven along the $`y`$-direction. By applying driving force $`\stackrel{}{𝒙}_\mathrm{C}`$ goes slowly to the position F, which is the point of inflection of the total potential. As soon as $`\stackrel{}{𝒙}_\mathrm{C}`$ arrives at F, $`\stackrel{}{𝒙}_\mathrm{C}`$ slips to the next sticking position close to G, which is the position of the total potential minimum. The sticking time around G is shorter than that around A due to elongation of the driving spring along the $`x`$-direction. The sticking times around A and G are about 70.0 and 30.0 % of the period of the stick-slip motion, respectively. As soon as $`\stackrel{}{𝒙}_\mathrm{C}`$ arrives at H, the flake slips again to the next sticking position close to J, which is equivalent to A. Then $`\stackrel{}{𝒙}_\mathrm{C}`$ repeats the periodic stick-slip motion and the trajectory of $`\stackrel{}{𝒙}_\mathrm{C}`$ becomes zigzag. The slipping times are of order of ten pico-seconds. The slips of $`\stackrel{}{𝒙}_\mathrm{C}`$ occur between positions close to the substrate potential minimum, in which the flake makes the AB stacking structure of the graphite, as shown in the inset of fig. 3. The atomic scale stick-slip motion results from the binding of the flake close to the AB stacking positions of the graphite flake on the graphite substrate.

Figure 4 shows the frictional force $`F`$ as a function of $`x_\mathrm{B}`$ for $`V_\mathrm{S}`$ $`=`$ 0.22 m/s driven along the $`x`$ (a) and $`y`$ (b) directions, respectively. $`F`$ initially increases linearly with $`x_\mathrm{B}`$ due to the elastic deformation of the driving spring. After $`F`$ takes the maximum value the curves exhibit periodic saw-tooth shapes, which are observed in experiments Mate ; Morita ; Hoshi ; Irr , due to the periodic stick-slip motions of the flake. The regions A, D and G as shown in the insets in fig. 4 (a,b) correspond to sticking states around those in fig. 3, respectively.

We define the maximum static frictional force $`F_s`$ and the kinetic frictional force $`F_k`$ as the maximum frictional force during driving and a time averaged value of the frictional force during one cycle of the stick-slip motion, respectively. Figure 5 shows $`F_s`$ and $`F_k`$ as a function of the loading force, $`L`$, driven along the $`x`$ and $`y`$-directions. The dependence of $`F_s`$ and $`F_k`$ on $`V_\mathrm{S}`$ is weak. This is because that even in the present calculation the time scale of the driving is much longer than that of atomic one. The values of the frictional force in fig. 5 are obtained by extrapolating data for $`V_\mathrm{S}`$ $`=`$ 0.22, 0.44 and 0.88 m/s to $`V_\mathrm{S}`$ $`=`$ 0 m/s for the comparison with experiments, in which the typical magnitude of $`V_\mathrm{S}`$ is about 40 nm/s. $`F_s`$ and $`F_k`$ linearly depend on the loading force and have a finite adhesive term $`a`$, which is the frictional force at $`L=0`$ nN. The frictional coefficient $`\mu `$ is defined as,

$`\stackrel{~}{F}=\mu L+a.`$ (3)
Here $`\stackrel{~}{F}`$ denotes $`F_s`$ or $`F_k`$. The magnitude of the frictional force in the range of the loading force in fig. 5 is very low with the frictional coefficients $`\mu _s`$ = $`\mu _k`$ $`=`$ 0.013 driven along the $`x`$-direction and $`\mu _s`$ $`=`$ $`\mu _k`$ $`=`$ 0.0055 driven along the $`y`$-directions. Here $`\mu _s`$ and $`\mu _k`$ denote the frictional coefficients of $`F_s`$ and $`F_k`$, respectively. These values are close to those observed in the experiment, $`\mu `$ = 0.012 Mate .

Here we discuss the mechanism of the low frictional coefficients of graphite. It has been believed that the low frictional coefficient of macroscopic graphite surfaces or graphite solid lubricants results from the weak interlayer interaction. We obtain, however, the frictional coefficient $`\mu `$ = 0.099 for a single carbon atom tip on the graphite substrate for the same pressure value with the case of the flake, which is much larger than those in our simulation of the flake and that in the experimentMate . The difference of the frictional coefficients indicates an another mechanism of the low frictional coefficients between graphite surfaces. In the flake there are two kinds of lattice sites $`\widehat{\mathrm{A}}`$ and $`\widehat{\mathrm{B}}`$ catching different forces from the substrate as shown in fig. 6.

The substrate potential for the flake $`V_{\mathrm{sub}}`$ is the sum of the substrate potential for a single atom at the sites $`\widehat{\mathrm{A}}`$, $`V_{\mathrm{sub}}^{\widehat{\mathrm{A}}}`$, and that at the sites $`\widehat{\mathrm{B}}`$, $`V_{\mathrm{sub}}^{\widehat{\mathrm{B}}}`$. That is, $`V_{\mathrm{sub}}`$ $`=`$ $`N_{\widehat{\mathrm{A}}}V_{\mathrm{sub}}^{\widehat{\mathrm{A}}}`$ $`+`$ $`N_{\widehat{\mathrm{B}}}V_{\mathrm{sub}}^{\widehat{\mathrm{B}}}`$, where $`N_{\widehat{\mathrm{A}}}`$ and $`N_{\widehat{\mathrm{B}}}`$ denote the number of the sites $`\widehat{\mathrm{A}}`$ and $`\widehat{\mathrm{B}}`$. In the flake $`N_{\widehat{\mathrm{A}}}`$ = $`N_{\widehat{\mathrm{B}}}`$. Here we consider the case that $`\stackrel{}{𝒙}_\mathrm{C}`$ moves straight between the nearest substrate potential minimums A and C shown in fig. 3. The atoms at sites $`\widehat{\mathrm{A}}`$ move from the potential minimum to the maximum and those at sites $`\widehat{\mathrm{B}}`$ move from the maximum to the minimum as shown in fig. 7.

During the motion the increase of $`V_{\mathrm{sub}}^{\widehat{\mathrm{A}}}`$ and the decrease of $`V_{\mathrm{sub}}^{\widehat{\mathrm{B}}}`$ cancel out each other. Then the substrate potential for the flake $`V_{\mathrm{sub}}`$ does not vary in the route between the nearest potential minimums. As a result the flat potential valleys appears for the flake in the shaded region shown in fig. 3. The valleys enable the flake to move easily and yield low frictional force. Due to this mechanism graphite becomes a good lubricant with low frictional coefficient. On the other hand the single carbon atom above the graphite substrate has large frictional coefficient because of the absence of such cancellation mechanism of the substrate potential.

The frictional coefficient mainly depends on the trajectory of $`\stackrel{}{𝒙}_\mathrm{C}`$. When the flake is large enough that the interaction between the flake and the substrate is so strong as to overcome the stiffness of the driving spring, the dominant contribution to the acceleration of the flake is flake size independent value, $`(_\mathrm{A}V_{\mathrm{sub}}^{\widehat{\mathrm{A}}}+_{\widehat{\mathrm{B}}}V_{\mathrm{sub}}^{\widehat{\mathrm{B}}})/2m_\mathrm{C}`$, for the same separation between the flake and the substrate. Here $`_{\widehat{\mathrm{A}},\widehat{\mathrm{B}}}`$ denote the coordinate derivatives at the sites $`\widehat{\mathrm{A}}`$ and $`\widehat{\mathrm{B}}`$, respectively. Then because the trajectory is independent of the flake size, the frictional coefficient is also independent of the flake size. The flake in this simulation is sufficiently large, so that the frictional coefficient is independent of the flake size. In fact the difference between the frictional coefficients for $`N_\mathrm{f}`$ = 42 and 110 is much small as about 2 %.

In the present work we have investigated atomic scale friction between clean graphite surfaces by numerical simulation. The simulation reproduces atomic scale stick-slip motion and low frictional coefficients of graphite. The former is due to the binding of the flake close to the stacking configurations on the substrate. The latter results from the cancellation of the forces coming from the substrate between two kinds of lattice sites in the flake, which is discovered in this letter for the first time within our knowledge.

We express sincere thanks to Profs. K. Miura and S. Morita for valuable discussions. The present work is financially supported by Grant-in-Aid for Scientific Research from Japan Society for the Promotion of Science.


text:
# The lack of variability of the iron line in MCG–6-30-15: general relativistic effects

## 1 Introduction

The X–ray spectrum of the Seyfert 1 galaxy MCG–6-30-15 exhibits a broad emission feature peaked at $`6.4`$ keV which extends from below $`4`$ keV to approximately $`7`$ keV, and was first resolved with ASCA \[Tanaka et al. 1995\]. The profile of this feature is consistent with that expected from iron fluorescence from the surface of an accretion disc surrounding a massive black hole.

The X–ray continuum emission of MCG–6-30-15 is highly variable (see e.g. Vaughan, Fabian & Nandra 2003; Fabian & Vaughan 2003). Since the reflection features, such as the observed iron line, are expected to be due to reprocessing of the hard X–ray continuum by dense gas in the accretion disc \[George & Fabian 1991, Matt, Perola & Piro 1991\], the line flux should respond to the changes in the illuminating continuum (Reynolds et al 1999). Line flux variations have been reported on timescales of about $`10^4`$ s, but with a much smaller amplitude than the observed continuum changes \[Iwasawa et al. 1996, Vaughan & Edelson 2001, Shih, Iwasawa & Fabian 2002, Fabian & Vaughan 2003\]. Furthermore, the variability of the iron line flux and the observed continuum are uncorrelated.

The source spectrum and variability can be explained by a two component model \[Shih, Iwasawa & Fabian 2002\] consisting of a variable power-law component and an almost constant reflection component containing the iron line. This model successfully accounts for the behaviour of MCG–6-30-15 during a $`320`$ ks long XMM–Newton observation \[Fabian & Vaughan 2003\]. Hereafter, we shall refer to these components as the Power Law Component (PLC) and the Reflection–Dominated Component (RDC). The existence of a varying soft component with a spectral shape uncorrelated with flux together with an almost constant harder component is also evident in the linearity of the flux–flux relationship recently presented by Taylor, Uttley & McHardy (2003).

The extreme red wing of the iron line indicates that the inner radius of the disc extends down to only a few gravitational radii (Wilms et al 2001; Fabian et al 2002). Strong gravitational field effects dominate the behaviour of photons in this region with not only gravitational redshift being severe but also gravitational light bending \[Martocchia & Matt 1996, Martocchia, Karas & Matt 2000, Dabrowski & Lasenby 2001, Martocchia, Matt & Karas 2002\]. Fabian & Vaughan (2003) have recently suggested that changes in the position of the power–law emission region could lead to large variations in the observed flux together with a strong and almost constant reflection component, containing the broad iron line. This behaviour would be due predominantly to strong light bending, which is expected if the primary source is located close to the central black hole and illuminates the inner regions of the accretion disc, as the shape of the iron line suggests. As the source approaches these regions of the disc close to the black hole so more radiation is bent away from our line of sight and intercepted by the disc. Such light bending also explains the high observed equivalent width of the line.

In this Letter, we investigate, for the parameters relevant to MCG–6-30-15, the variability of the PLC and of the RDC as a function of the position and state of motion of the primary X–ray source, seeking solutions that allow for a strongly varying PLC and an almost constant RDC. The main purpose of this work is to provide an astrophysical context in which the two component model for the variability of MCG–6-30-15 can be explained self–consistently. Our calculations are performed using a combination of Monte Carlo and ray-tracing methods that take into account the effects of Doppler/gravitational redshifts and light beaming/bending on the motion of photons in the spacetime of a rotating black hole.

## 2 Model and assumptions

The observed broad iron line profile and its extended red wing \[Wilms et al. 2001, Fabian et al. 2002\] imply that most of the emission originates at radii less than $`6r_g`$ ($`r_g=GM/c^2`$), suggesting that the central black hole is rapidly spinning \[Iwasawa et al. 1996, Dabrowski et al. 1997\]. Our computations are then performed in the spacetime of a maximally rotating (Kerr) black hole.

The accretion disc is assumed to be thin and to lie in the equatorial plane, perpendicular to the rotation axis. The accreting material is flowing along stable circular geodesics and the disc extends from the radius $`r_{\mathrm{ms}}`$ of the marginal stable circular orbit \[e.g. Novikov & Thorne 1973\] out to $`r_{\mathrm{out}}=100r_g`$.

The source of primary hard X–rays is assumed to have a ring–like axisymmetric geometry and is located above the disc at a distance $`\rho _s`$ from the rotation axis and at height $`h_s`$, related to the Boyer–Lindquist coordinates $`r_s`$ and $`\theta _s`$ of the source by $`r_s=\sqrt{\rho _s^2+h_s^2}`$ and $`\theta _s=\mathrm{arctan}(\rho _s/h_s)`$. The source can be both static or corotating with the disc. If the source is static, its $`4`$–velocity has to be proportional to the time–like Killing vector of the spacetime so that

$$𝐮=C_{\mathrm{stat}}_t,$$
(1)
while if it is corotating one has

$$𝐮=C_{\mathrm{corot}}(_t+\mathrm{\Omega }_\varphi ),$$
(2)
where $`\mathrm{\Omega }1/[a+(r\mathrm{sin}\theta )^{3/2}]`$ and where $`C_{\mathrm{stat}}`$ and $`C_{\mathrm{corot}}`$ are found by requiring that the source follows a time–like world line ($`𝐮𝐮=1`$).

The case of a static source is investigated mainly for comparison with the corotating one that we consider physically more plausible. This is because since flares are believed to be associated with magnetic activity of the disc, they are more likely to be corotating with the accreting material rather than static. In the case of a corotating source, since the orbital timescale in the vicinity of the central massive black hole is much shorter than the integration time needed in observations, any information on the azimuthal position of the flare is lost. This justifies our choice of a ring–like axisymmetric source rather than a point–like whose azimuthal position would affect the iron line profile emitted from the disc but is unobservable with current long observations.

The source produces isotropic emission in its proper frame with a power law luminosity $`L_s=L_0E_s^\alpha `$, where $`E_s`$ is the photon energy in the source comoving frame and the spectral index is chosen to be $`\alpha =1.1`$, consistent with the observed photon index in the spectrum of MCG–6-30-15 (see e.g. Fabian & Vaughan 2003). Isotropic emission is enforced by using the HEALPix111see http://www.eso.org/science/healpix (Hierarchical Equal Area isoLatitude Pixelisation) package which allows one to define equal area curvilinear quadrilaterals on the sphere. Isotropy (in the source proper frame) is obtained by sending photons along the directions defined by the geometrical centre of each quadrilateral.

The photons emitted by the primary source illuminate both the accretion disc and the collecting area of a distant observer, whose distance from the source is taken to be $`r_{\mathrm{obs}}=10^3r_g`$ for numerical purposes. The observer measures both the direct flux and the reflected emission from the disc. The observer inclination has been fixed to $`30`$ degrees, which is appropriate for the case of MCG–6-30-15.

In this work, we assume that the direct continuum emission from the primary source dominates over the reflected continuum emission, which is neglected in our computations. However, since the reflected continuum varies together with the iron line, it does not affect the results on the variability of the PLC and of the RDC (hereafter represented by the iron line only), which are the main purpose of this Letter.

Each simulation is carried out by sending a given number $`N`$ of photons from the primary source and by integrating along the null geodesics in the Kerr spacetime until the photons reach the accretion disc (or are lost into the black hole horizon) or the observer collecting area. This procedure allows us to compute directly the PLC at the observer and the illuminating flux on the disc as a function of the photons incident angle $`\theta _i`$ and energy $`E_d`$, taking into account both special and general relativistic effects.

The iron line local intensity is computed by assuming that the fluorescent photons are produced by cold, non–ionised matter at the rest frame energy of $`6.4`$ keV and by making use of the work by George & Fabian (1991) who provided analytical approximations for the dependence of the fluorescent emission on $`\theta _i`$ and $`E_d`$ (see also Ruszkowski 2000; Lu & Yu 2001 for further details). In this work, the iron line emission is assumed to be isotropic in the rest frame of the disc, which is not a strong limitation for low inclination objects such as MCG–6-30-15.

To calculate line profiles and fluxes, we use a ray–tracing technique to follow the trajectories of photons from the observer until they either intersect the accretion disc or the black hole horizon and we calculate the redshift factor corresponding to each arrival position on the disc \[Laor 1991, Dabrowski et al. 1997\]. Then the line flux is obtained by making use of the relativistic invariance of $`I_\nu /\nu ^3`$, where $`I_\nu `$ is the specific line intensity, and by integrating over the solid angle subtended by the whole disc at the observer.

## 3 PLC and RDC variability

We investigate the effects that the position and the motion of the primary source have on the continuum and reflected observed fluxes, with particular emphasis on the associated variability of both the PLC and RDC components.

As already pointed out in previous studies \[Martocchia & Matt 1996, Martocchia, Karas & Matt 2000, Martocchia, Matt & Karas 2002, Dabrowski & Lasenby 2001\], if the primary source is very close to the central massive black hole, a large fraction of the emitted photons will be bent onto the disc (or lost into the hole event horizon), reducing the observed direct flux with respect to the reflected one and enhancing the iron line equivalent width (EW).

This effect is shown in Fig. 1 where we consider the arrival positions of the photons emitted from a static source located at $`\rho _s=2.0r_g`$ as a function of the height ($`h_s`$) of the source. As an example, for a source height of $`h_s=2r_g`$, the percentage of photons that reach the accretion disc is 64 per cent. The remaining 36 per cent is equally distributed between photons that are lost in the hole horizon and those that reach the observer at infinity. Similar results are obtained also in the case of a corotating source. However in the latter case, fewer photons are lost into the hole horizon (8 per cent) due to the strong beaming in the direction of the source motion that reduces the number of photons directed towards the hole. The illumination pattern on the disc is then more extended for a corotating source (74 per cent of the photons hit the disc) than for a static one.

Here, the PLC variability is supposed to be caused by variations of the primary source position. Since the source height controls the amount of direct radiation that reaches the observer (see Fig. 1), the observed PLC from a source close to the black hole is much smaller than that of a source located at larger height. At the same time, the RDC responds to a much larger region than the PLC, especially if the source is corotating, and it is likely to vary less than the power law component.

To illustrate this behaviour, and to compare the variation of the two components we present results for sources located at $`\rho _s=2r_g`$ from the rotation axis, but similar results can be obtained within $`\rho _s34r_g`$. In Figs. 2 and 3 we show the PLC and RDC fluxes as a function of the source height for a static and a corotating source respectively. The PLC flux is defined as the observed flux of the power law component evaluated at $`6.4`$ keV, while the RDC is the total iron line observed flux (i.e. the integral over energy of the line). We also show the iron line equivalent width defined as the ratio between the RDC and the PLC. The general behaviour is that the variation in EW is dominated by the PLC variability and is clearly anti–correlated with flux. Since we neglect the contribution of the reflected continuum, the values of the line EW are overestimated, especially in the most extreme cases (low source height) and they have to be taken only as an indication.

Both the power law and the reflection–dominated components vary with the height of the primary source, but the variation of the RDC has a much smaller amplitude. This behaviour is enhanced in the corotating case shown in Fig. 3 where, if the source height changes from $`3`$ to $`8r_g`$, a variation of the PLC by almost a factor $`4`$ is accompanied by a variation in the RDC with a maximum amplitude of 15 per cent around its mean value. Even larger variations of the PLC are found for smaller values of $`\rho _s`$, together with an almost constant RDC. These results can be compared with the work by Fabian & Vaughan (2003) who analysed the spectral variability of MCG–6-30-15 finding variations of the RDC up to about 25 per cent within a PLC change of about a factor $`4`$.

We note here that the two component model also accounts (Vaughan et al 2003, in preparation) for the shape and amplitude of the r.m.s. variability spectrum of MCG–6-30-15, plotted in Fabian et al. (2002).

### 3.1 Emissivity and line profiles

We compute the local iron line emissivity profile induced by a ring–like primary source by rotating the illumination pattern of the one for a point–like source over the azimuthal position on the disc. The emissivity is thus axisymmetric and it is a function of the radial position on the disc only. The emissivity $`ϵ(r)`$ produced by a corotating ring source at $`\rho _s=2r_g`$ is shown in Fig. 4 for four different source heights. Due to the anisotropy of incident radiation (controlled by $`h_s`$), the emissivity in the inner regions is reduced by increasing the source height, while it increases in the outer disc (see also Martocchia, Karas & Matt 2000).

At a fixed $`h_s`$, the emissivity is steeper in the inner disc than in the outer regions with a transition region between about $`3`$ and $`10r_g`$. This is in good agreement with a recent analysis \[Fabian et al. 2002\] where a best fit to the iron line of MCG–6-30-15 was found by using a broken power law emissivity with a break radius of about $`6r_g`$. The outer emissivity index was found to be $`2.5`$, while a steeper index of about $`4.8`$ was derived within $`6r_g`$.

The iron line profiles are broad with an extended red wing if the primary source is located close to the black hole, while the profile becomes progressively narrower as the distance from the hole is increased. This is due to the different illumination pattern on the accretion disc (and different emissivity) for different source heights. If the source is low, the inner regions of the disc are more illuminated resulting in the enhancement of the relativistic broadening of the iron line. Our model then predicts broader line profiles when the source flux is low (source at low heights). This is indeed the sense of the changes seen in the ASCA 1994 observation \[Iwasawa et al. 1996\]. In Fig. 5 we show the line profiles obtained in the case of a corotating source at $`\rho _s=2.0r_g`$ for different heights in the range $`38r_g`$, whose variation accounts for a change in the PLC by a factor $`4`$.

The changes in the line profiles for different source heights are subtle and it could be difficult to detect them in real X-ray data. If this is the case, an observer would conclude that the line is almost constant both in flux and in spectral shape while the PLC varies in normalisation, as it seems to be the case in the recent XMM-Newton observation of MCG–6-30-15 \[Fabian et al. 2002\] that was used to study the source variability \[Fabian & Vaughan 2003\]. This possibility has been tested by comparing the line profiles we produced with this same data set. We analysed the $`3`$$`10`$ keV time–averaged spectrum obtained with the EPIC pn camera using XSPEC version 11 \[Arnaud 1996\]. The spectrum was fitted using a simple power law plus the relativistic line produced by our model when a corotating source is placed at $`\rho _s=2.0r_g`$ with variable height between $`3`$ and $`8r_g`$.

The $`\chi ^2`$ difference between the two more extreme heights is $`\mathrm{\Delta }\chi ^213`$ with $`\chi ^2=265.8`$ for $`h_s=3r_g`$ and $`\chi ^2=252.9`$ for $`h_s=8r_g`$, both with $`178`$ degrees of freedom. However, this simple model leaves some residuals around $`6.4`$ keV and the significance of the preferred source height is poor. These residuals can be accounted for by an additional narrow Gaussian component with an equivalent width of only $`30`$ eV that can be due to reflection from distant material and is not included in our model for line emission from the inner $`100r_g`$ of the accretion disc. With the addition of this component, the emission line in the XMM-Newton data is well fitted by our model and, most remarkably, the $`\chi ^2`$ for the most extreme heights is now identical (both models have $`\chi ^2/\nu =190.5/1751.09`$).

This means that, once the narrow line component is included, the changes in the line shape are substantially reduced and may be difficult to detect with current instruments. If the narrow component is emitted from distant material, it is likely to remain constant. Future observations of the source may be capable of resolving the narrow component from the broader one emitted from the accretion disc, allowing more stringent constraints to be placed on the height of the illuminating primary source. In summary, the model we have explored consistently accounts for the uncorrelated variability of the PLC and RDC components and produces also iron line profiles in good agreement with the most recent X–ray observation of MCG–6-30-15.

## 4 Conclusions

The spectral variability of MCG–6-30-15 can be explained by an almost constant Reflection–Dominated Component together with a highly variable Power Law Component.

Here we have explored a model in which the variability of the PLC is accounted for by changes in the height of the primary source above the accretion disc. Light bending by the gravitational field of the black hole substantially reduces the PLC flux at low heights, allowing for high continuum variability. The iron line flux, which represents the RDC, is much more constant with an amplitude of variation of about 15–20 per cent when the PLC varies by a factor $`4`$.

It should be noted that the region where the flare (or multiple flares) creating the PLC should be located to reproduce this behaviour is close to the rotation axis, within $`3`$-$`4`$ gravitational radii. It is most plausibly powered by magnetic fields from the disc, possibly linking to the hole (e.g. Blandford & Znajek 1977; Wilms et al 2001), and perhaps forming the base of a jet (note that the radio emission from MCG–6-30-15 is weak and unresolved, Morganti et al. 1996).

We also find that the iron line emissivity on the accretion disc is consistent with a broken power law profile with a transition between about $`3`$ and $`10r_g`$, in good agreement with recent observations.

In summary, the strong light bending close to the black hole can cause the PLC and RDC to appear disconnected. Even if the source of the PLC is intrinsically constant in luminosity, it will appear to an outside observer to vary if its height from the disc changes. Much of the radiation is bent onto the disc in a manner which means that the flux of the RDC will appear constant. Subtle profile changes are expected in the sense that the line will appear narrower when the PLC is bright and broader when dim, but they are difficult to detect with current instruments.

It is of course possible that the source of the PLC undergoes intrinsic luminosity variations on both short and long timescales. These will lead to connected variability between the PLC and RDC, and reverberation on short timescales. For MCG–6-30-15 these timescales are probably less than $`10^3`$ s, too short to be investigated with XMM-Newton but within the grasp of Constellation-X.

## Acknowledgements

We would like to thank the referee, Dr. Giorgio Matt, for his helpful comments. Some of the results have been derived using the HEALPix package \[Górski, Hivon & Wandelt 1999\]. ACF thanks the Royal Society for support.


text:
# Energy spectra of the ocean’s internal wave field: theory and observations.

## Abstract

The high-frequency limit of the Garrett and Munk spectrum of internal waves in the ocean and the observed deviations from it are shown to form a pattern consistent with the predictions of wave turbulence theory. In particular, the high frequency limit of the Garrett and Munk spectrum constitutes an exact steady state solution of the corresponding kinetic equation.

-

Introduction.  Internal waves are an important piece of energy and momentum budgets for the earth’s atmosphere and ocean. The drag associated with internal wave breaking needs to be included in order to obtain accurate simulations of the atmospheric Jet Stream  and it has been argued that the ocean’s Meridional Overturning Circulation  is forced by the diffusion of mass  associated with internal wave breaking  rather than by the production of cold, dense water by convection at high latitudes. Both circulations represent important pieces of the earth’s climate system.

In a classical work , Garrett and Munk demonstrated how observations from various sensor types could be synthesized into a combined wavenumber-frequency spectrum, now called the Garrett-and-Munk (GM) spectrum of internal waves. Consistent only with linear internal wave kinematics, the GM spectrum was developed as an empirical curve fit to available data. Even though deviations have been noted near boundaries , and at the equator , the last significant model revision  has surprisingly stood the test of time. However, a review of open ocean data sets reveals subtle variability in spectral power laws. We show in this letter that predictions based upon a weakly nonlinear wave turbulence theory are consistent with both the high frequency–high wave number limit of GM spectrum and the observed variability.

In this letter, we will consider only the high frequency–high wave number limit of GM; for brevity, we shall denote this henceforth as GMh. The GMh is given by

$$E(m,\omega )Nm^2\omega ^2.$$
(1)
Here $`E`$ is the spectral wave energy density, $`N`$ the buoyancy frequency, $`m`$ the vertical wavenumber, and $`\omega `$ the frequency. The total energy density of a wavefield is $`E=E(m,\omega )𝑑m𝑑\omega `$.

The possibility that the internal wavefield might exhibit a universal character represents an attractive theoretical target, and much effort (as reviewed by ) was devoted to studying the issue of nonlinearity in the context of resonant wave interactions. That line of work is based on a Lagrangian description of the flow, with two main approximations: that fluid particles undergo small displacements, and that nonlinear interactions take place on a much longer time–scale than the underlying linear dynamics. An approximate kinetic equation describing the time evolution of spectral wave energy was derived, and it was shown  that the GMh spectrum (1) was close to being a stationary solution.

An alternative to the Lagrangian formulation, based on a Hamiltonian description in isopycnal (density) coordinates, was recently proposed . This approach does not invoke a small–displacement assumption and yields a comparatively simple kinetic equation with an exact steady power–law solution in the high frequency limit. That steady state solution, \[see (12) below\] is close to the GMh spectrum (1), yet there is a noticeable difference. Motivated by this difference, we tried to estimate the accuracy of the GMh power laws and thus reviewed extant observations from the literature. In the process of analyzing the data, we found that there was subtle variability in the high wavenumber, high frequency spectrum, forming a distinct pattern.

We then reexamined the kinetic equation of  and found its full family of steady state solutions, of which the solution reported in  is just one member. This family of solutions compares well with the variability found in the observations. Moreover, the GMh spectrum (1) is a member of this family, thus describing the GMh spectrum simply as an exact steady–state solution to the kinetic equation derived in .

Hence, in this article we present evidence for variability in the high frequency-high wavenumber open-ocean internal wavefield, and find that a wave turbulence approach predicts that both GMh itself and the observed variability are stationary states of the kinetic equation. The variability itself, and its likely roots in variable forcing, Coriolis effects, underlying stratification and currents, as well as the low frequency range of the energy spectrum, are fundamental problems posing exciting challenges for future research.

Overview of observations: a family of spectra.  Below we present a summary of historical oceanic internal wave energy spectra. These observations are re-analyzed to study whether the high-frequency, high-wavenumber spectra may form a pattern. We review seven data sets available in the literature. We shall present a detailed analysis of these data sets elsewhere; here we just list them along with their high–frequency, high–wave number asymptotics. Let us assume that, in this limit, the three dimensional wave action $`n(𝐤,m)`$ can be approximated by horizontally isotropic power laws of the form

$$n_{𝐤,m}=n_0|𝐤|^x|m|^y,$$
(2)
where $`𝐤`$ is the horizontal wave vector, $`k=|𝐤|`$ its modulus, $`m`$ vertical wavenumber and $`n_0`$ is a constant.

Using the linear dispersion relation of internal waves, $`\omega _{𝐤,m}|𝐤|/m`$, this action spectrum can be transformed from the wavenumber space $`(k,m)`$ to the vertical wave-number-frequency space $`(\omega ,m)`$ . Multiplication by the frequency yields the corresponding energy spectrum,

$$E(m,\omega )\omega ^{2x}m^{2xy}.$$

The total energy of the wave field is then

$$E=\omega (𝐤,m)n(𝐤,m)𝑑𝐤𝑑m=E(\omega ,m)𝑑\omega 𝑑m.$$

Below we list extant data sets with concurrent vertical profile and current meter observations and some major experiments utilizing moored arrays, along with our best estimate of their high-wave-number high frequency asymptotics (the order is chronological):
$``$ The Mid-Ocean Dynamics Experiment (MODE), March-July 1973, Sargasso Sea ($`26^{}0^{}`$ N, $`69^{}40^{}`$ W): $`m^{2.25}\omega ^{1.6}`$ ;
$``$ The Internal Wave Experiment (IWEX), 40 days observations in November-December 1973, Sargasso Sea thermocline ($`27^{}44^{}`$ N, $`69^{}51^{}`$ W): $`k^{2.4\pm 0.4}\omega ^{1.75}`$ ;
$``$ The Arctic Internal Wave Experiment (AIWEX), March to May of 1985, Canada Basin thermocline, ($`74^{}`$ N, $`143146^{}`$ W): $`m^{2.15}\omega ^{1.2}`$ ;
$``$ The Frontal Air-Sea Interaction Experiment (FASINEX), January to June of 1986, Sargasso Sea thermocline ($`27^{}`$ N, $`70^{}`$ W): $`m^{1.9\mathrm{to}2.0}\omega ^{1.75}`$ ;
$``$ Patches Experiment (PATCHEX), 7.5 days during October of 1986, eastern Subtropical North Pacific, ($`34^{}`$ N, $`127^{}`$ W): $`m^{1.75}\omega ^{1.65\mathrm{to}2.0}`$ ;
$``$ The Surface Wave Process Program (SWAPP) experiment, 12 days during March, 1990, eastern Subtropical North Pacific thermocline, ($`35^{}`$ N, $`127^{}`$ W): $`m^{1.9}\omega ^{2.0}`$ ;
$``$ North Atlantic Tracer Release Experiment (NATRE), February-October 1992, eastern Subtropical North Atlantic thermocline, ($`26^{}`$ N $`29^{}`$ W): $`m^{2.75}\omega ^{0.6}`$ (for $`1<\omega <6`$cpd) .
These deep ocean observations (Figure 1) exhibit a higher degree of variability than one might anticipate for a universal spectrum. Moreover, the deviations from the GMh spectral power laws form a pattern: they seem to roughly fall upon a curve with negative slope in the $`(x,y)`$ plane. We show in the next section that the predictions of wave turbulence theory are consistent with this pattern.

A wave turbulence formulation for the internal wave field. In this section we assume that the internal wave field can be viewed as a field of weakly interacting waves, thus falling into the class of systems describable by wave turbulence. Wave turbulence is a universal statistical theory for the description of an ensemble of weakly interacting particles, or waves. This theory has contributed to our understanding of spectral energy transfer in complex systems , and has been used for describing surface water waves since pioneering works by Hasselmann , Benney and Newell  and Zakharov .

The dynamics of oceanic internal waves can be most easily described in isopycnal (i.e. density) coordinates, which allow for a simple and intuitive Hamiltonian description . To describe the wave field, we introduce two variables: a velocity potential $`\varphi (𝐫,\rho )`$, and an isopycnal straining $`\mathrm{\Pi }(𝐫,\rho )`$. The horizontal velocity is given by the isopycnal gradient $``$ of the velocity potential, $`u(𝐫,\rho )=\varphi (𝐫,\rho )`$. The straining $`\mathrm{\Pi }=\rho /(_z\rho )`$ can also be interpreted as the fluid density in isopycnal coordinates.

These two variables form a canonically conjugated Hamiltonian pair, so that the primitive equations of motion (i.e. conservation of horizontal momentum, hydrostatic balance, mass conservation and the incompressibility constraint) can be written as a pair of canonical Hamilton’s equations,

$`_t\mathrm{\Pi }=\delta /\delta \varphi ,_t\varphi =\delta /\delta \mathrm{\Pi },`$ (3)
$`={\displaystyle \frac{1}{2}}{\displaystyle \left(\mathrm{\Pi }|\varphi |^2\left|^\rho \frac{\mathrm{\Pi }}{\rho _1}𝑑\rho _1\right|^2\right)𝑑𝐫𝑑\rho }.`$ (4)
(5)
The first term in the Hamiltonian clearly corresponds to the kinetic energy of the flow, the second term can be shown to correspond analogously to the potential energy.

Performing the Fourier transform, and introducing a complex field variable $`a_𝐩`$ via

$$\varphi _𝐩=\frac{iN\sqrt{\omega _𝐩}}{\sqrt{2g}k}\left(a_𝐩a_𝐩^{}\right),\mathrm{\Pi }_𝐩=\frac{\sqrt{g}k}{\sqrt{2\omega _𝐩}N}\left(a_𝐩+a_𝐩^{}\right),$$

where $`𝐩=(𝐤,m)`$ is the 3-d wave vector, $`N`$ is buoyancy frequency, $`g`$ is gravity acceleration; the canonical pair of equations of motion and the Hamiltonian (5) read

$`i{\displaystyle \frac{}{t}}a_𝐩={\displaystyle \frac{}{a_𝐩^{}}},\text{where}={\displaystyle \omega _𝐩|a_𝐩|^2𝑑𝐩}+`$ (6)
$`{\displaystyle V_{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}}^{𝐩_\mathrm{𝟑}}\left(a_{𝐩_\mathrm{𝟏}}^{}a_{𝐩_\mathrm{𝟐}}^{}a_{𝐩_\mathrm{𝟑}}+a_{𝐩_\mathrm{𝟏}}a_{𝐩_\mathrm{𝟐}}^{}a_{𝐩_\mathrm{𝟑}}^{}\right)\delta _{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}𝐩_\mathrm{𝟑}}𝑑𝐩_\mathrm{𝟏}𝑑𝐩_\mathrm{𝟐}𝑑𝐩_\mathrm{𝟑}}+`$ (7)
$`{\displaystyle \frac{1}{6}V_{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}}^{𝐩_\mathrm{𝟑}}\left(a_{𝐩_\mathrm{𝟏}}^{}a_{𝐩_\mathrm{𝟐}}^{}a_{𝐩_\mathrm{𝟑}}^{}+a_{𝐩_\mathrm{𝟏}}a_{𝐩_\mathrm{𝟐}}a_{𝐩_\mathrm{𝟑}}\right)\delta _{𝐩_\mathrm{𝟏}+𝐩_\mathrm{𝟐}+𝐩_\mathrm{𝟑}}𝑑𝐩_\mathrm{𝟏}𝑑𝐩_\mathrm{𝟐}𝑑𝐩_\mathrm{𝟑}},`$ (8)
with wave-wave interaction matrix elements given by : $`V_{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}}^𝐩=U_{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}}^𝐩+U_{\mathrm{𝐩𝐩}_\mathrm{𝟐}}^{𝐩_\mathrm{𝟏}}+U_{\mathrm{𝐩𝐩}_\mathrm{𝟏}}^{𝐩_\mathrm{𝟐}}`$ with

$$U_{𝐩_\mathrm{𝟏}𝐩_\mathrm{𝟐}}^𝐩=\frac{N}{4\sqrt{2g}}\frac{𝐤_\mathrm{𝟐}𝐤_\mathrm{𝟑}}{k_2k_3}\sqrt{\frac{\omega _{𝐩_\mathrm{𝟐}}\omega _{𝐩_\mathrm{𝟑}}}{\omega _{𝐩_\mathrm{𝟏}}}}k_1.$$

These field equations are equivalent to the primitive equations of motion for internal waves (up to the hydrostatic balance and Boussinesq approximation); the work reviewed in  instead resorted to a small displacement approximation to arrive at similar equations. We will argue elsewhere that this extra assumption does not provide an internally consistent description of interactions between extremely scale separated waves. For the purposes of this letter, it suffices to note that the two kinetic equations are different and yield different steady solutions.

We shall characterize the field of interacting internal waves by its wave action $`\delta _{𝐩𝐩^{}}n_𝐩=a_𝐩a_𝐩^{}^{}.`$

Under the assumption of weak nonlinear interaction, one derives a closed equation for the evolution of the wave action, the kinetic equation. Assuming horizontal isotropy, the kinetic equation can be reduced further by averaging over all horizontal angles, obtaining \[with $`p=(k,m)`$ and $`dp_1dp_2=dk_1dm_1dk_2dm_2`$\]

$`{\displaystyle \frac{dn_{k,m}}{dt}}={\displaystyle \frac{1}{k}}{\displaystyle \left(R_{p_1p_2}^pR_{pp_2}^{p_1}R_{p_1p}^{p_2}\right)𝑑p_1𝑑p_2/\mathrm{\Delta }_{k_1k_2}^k},`$ (9)
$`R_{p_1p_2}^p=\delta _{\omega _p\omega _{p_1}\omega _{p_2}}f_{p_1p_2}^p|V_{p_1p_2}^p|^2\delta _{mm_1m_2}kk_1k_2,`$ (10)
(11)
where $`f_{p_1p_2}^p=n_{p_1}n_{p_2}n_p(n_{p_1}+n_{p_2})`$ and $`\mathrm{\Delta }_{k_1k_2}^k=\left(2\left[(kk_1)^2+(kk_2)^2+(k_1k_2)^2\right]k^4k_1^4k_2^4\right)^{1/2}/2`$.

A family of steady state power-law solutions to the kinetic equation. In wave turbulence theory, three-wave kinetic equations admit two classes of exact stationary solutions: thermodynamic equilibrium and Kolmogorov flux solutions, with the latter corresponding to a direct cascade of energy –or other conserved quantities– toward the higher modes. The fact that the thermodynamic equilibrium –or equipartition of energy– $`n_𝐩=1/\omega _𝐩`$ is a stationary solution of (11) can be seen by inspection, whereas in order to find Kolmogorov spectra one needs to be more elaborate. In  we used the Zakharov-Kuznetsov conformal mapping  to show analytically that the following wave action spectrum constitutes an exact steady state solution of (11) \[note the difference with (1)\]:

$`n_{𝐤,m}=n_0|𝐤|^{7/2}|m|^{1/2};E(m,\omega )\omega ^{1.5}m^2,`$ (12)
Remarkably though, this is not the only steady state solution of the kinetic equation having nonzero spectral energy fluxes. In fact, there is a full family of such power law steady solutions. To see this, consider the kinetic equation (11), and substitute into it the ansatz (2). Let us now denote the resulting RHS of (11) by $`I(k,m)`$. For steady states, $`I(k,m)`$ needs to vanish for all values of $`k`$ and $`m`$, for appropriately chosen values of $`(x,y)`$. However, once $`I`$ vanishes for one such wavenumber $`(k,m)`$, it does so for all, due to the fact that $`I`$ is a bi–homogeneous function of $`k`$ and $`m`$:

$$I(\alpha k,\beta m)=\alpha ^{4+2x}\beta ^{1+2y}I(k,m).$$
(13)
Hence we can fix $`k`$ and $`m`$, and seek zeros of $`I`$ as a function of $`x`$ and $`y`$. The exact analytical solution (12) cannot correspond to an isolated zero of $`I`$, since $`(_xI,_yI)`$ is nonzero (it is proportional to the energy flux in the Kolmogorov solution ). Hence, by the Implicit Function Theorem, there must exist a curve of zeros of $`I(x,y)`$.

Since this family of steady state solutions is not all apparently amenable to a closed form, we sought the zeros of $`I`$ by numerical integration. This involves a certain amount of work. First, the delta-functions in (11) restrict contributions to the resonant set. Consider, for example, the resonant set

$$𝐤=𝐤_\mathrm{𝟏}+𝐤_\mathrm{𝟐},m=m_1+m_2,\omega _{k,m}=\omega _{k_1,m_1}+\omega _{k_2,m_2}.$$

Given $`𝐤`$, $`𝐤_\mathrm{𝟏}`$, $`𝐤_\mathrm{𝟐}`$ and $`m`$, one can find $`m_1`$ and $`m_2`$ satisfying this resonant condition by solving

$$k/m=k_1/|m_1|+k_2/|mm_1|.$$

This equation reduces to a quadratic equation for $`m_1`$, and then one can find $`m_2`$ from $`m_2=mm_1`$. After this reduction, one is left with a two–dimensional integral, over $`|k_1|`$ and $`|k_2|`$. This infinite domain is further restricted by the requirement that $`|k_1|`$, $`|k_2|`$ and $`|k|`$ are such that they can correspond to the sides of a triangle; this restricted (though still infinite) domain is called the *kinematic box* in the oceanographic literature. The next problem for the numerical integration is that the integrand diverges ( typically in an integrable fashion) at the boundaries of the domain. This is solved by a suitable change of variables. Finally, a second substitution renders the domain of integration finite.

The resulting family of zeros is depicted in Figure 1. Notice that the curve passes through the exact solution (12). More importantly, it also passes through the point $`(4,0)`$, corresponding to the GMh spectrum (1). Hence this classical spectrum is for the first time shown to correspond to an exact steady solution to a kinetic equation based on first fluid principles.

Finally, we note the integrals converge in the parameter regime occupied by the observations. In regions of tightly spaced contour lines ($`x<1.7`$ and $`y<0.7`$, $`x>4.2`$ and $`y<0.4`$) (4) is nonintegrable.

The other points marked on the figure correspond to the observational sets discussed above. Notice that, with the exception of NATRE, they all lie very close to the zeros of $`I`$. Therefore the predictions of wave turbulence are consistent with the observed deviations from GMh.

In fact, the NATRE point lies in an area of $`(x,y)`$ space where $`z=I`$ and $`z=0`$ are nearly tangential, thus making the line of zeros effectively “thicker” (in other words, the collision integral is not zero at the observed points, but it is very small, possibly allowing other, typically smaller effects to take over.)

– Conclusions We have shown that the wave turbulence formalism captures much of the variability apparent in the oceanic internal wave field. This includes the characterization of the spectral curve put together by Garrett and Munk as an exact steady solution to a kinetic equation for the evolution of the wave field, derived from first principles. In addition, the curve of steady solutions to this kinetic equation is consistent with much of the observed variability in the energy spectra. We conjecture that the placement along this curve of individual observations depends on the nature of the forcing (for instance, by tides and atmospheric winds), the local degree of stratification, vorticity and shear, and the variable magnitude of the Coriolis parameter. This is the subject of ongoing research.

Acknowledgments YL is supported by NSF CAREER grant DMS 0134955 and by ONR YIP grant N000140210528; KP is supported by NSF grant OCE 9906731; ET is supported by NSF grant DMS 0306729.


text:
# A novel view of plane wave expansion method in photonic crystals

## Abstract

We propose a method derived from the simple plane wave expansion that can easily solve the interface problem between vacuum and a semi-infinite photonic crystal. The method is designed to find the complete set of all the eigenfunctions, propagating or evanescent, of the translation operators $`\{𝐓_𝐑\}`$, at a fixed frequency. With these eigenfunctions and their eigenvalues, the transmitted and reflected waves can be determined. Two kinds of applications are presented for 2D photonic crystals. The first is a selection rule for determine the normal direction of the vacuum-photonic crystal interface to achieve the highest attenuation effect at a gap frequency. The second is to calculate the transmittance and reflectance for a light incident from vacuum to an semi-infinite photonic crystal. As an example we recalculate a system studied previously by K. Sakoda et al. and get results in agreement with theirs.

photonic crystal, plane wave, evanescent wave, transmittance, penetration length, numerical simulation  preprint: APS/123-QED

Since Yablonovitch 1 ; 2 ; 3  discovered a periodic dielectric structure that has an absolute gap in the frequency spectrum for electromagnetic waves, the idea of photonic crystals have attracted great interest. Many phenomena have been predicted theoretically and many application possibilities have been explored 2 ; 4 ; 5 ; 6 . Corresponding studies in the early years most authors paid their attention on the frequency spectrum gaps, and the most convenient method to calculate the band gaps is the plane wave expansion method Sakoda ; Gu1 . Recently, various kinds new methods have been proposed to compute some other relevant physical parameters such as transmittance and penetration depth 7 ; 8 ; 9 ; 17  for a finite system.

In this paper we also address the transmittance and penetration depth problems, but use a different method that includes all the information of propagating and evanescent modes getting from the translation operator. We show that by appropriately modifying the conventional plane wave expansion method we can enlarge its application region, and which makes it easy to solve the interface problem between vacuum and a semi-infinite photonic crystal system.

Our method has several advantages. First, both the “air rods in dielectric” and “dielectric rods in air” problems can be solved, without any restriction on the shapes of the rods and position of cutting plane that separating the semi-infinite photonic crystal region from the air region, which is impossible for the LKKR method 12 ; 13 ; 14 ; 15 ; 16 . Second, all information getting from the complete set of the eigenfunctions of the translation operator are used, including both the propagating and evanescent modes. This makes it easy to analyze and discuss phenomena using the well established knowledge of solid state physics. Third, the finite size effects such as the resonance behavior of the transmittance curve caused by the finite thickness of the photonic crystal sample can be easily isolated. We can thus accurately calculate the transmittance and reflectance for a very thick photonic crystal sample.

For a system without free charge and current, and if the permittivity $`ϵ(𝐫)`$ and permeability $`\mu (𝐫)`$ are scalars independent of time, the magnetic field $`𝐇(𝐫,t)`$ satisfies

$$\times \frac{1}{ϵ}\times 𝐇(𝐫,t)=\mu \frac{^2}{t^2}𝐇(𝐫,t),$$
(1)
where $`𝐇(𝐫,t)=\underset{\omega }{}𝐇_\omega (𝐫)e^{i\omega t}`$, and

$$\times \frac{1}{ϵ}\times 𝐇_\omega (𝐫)=\mu \omega ^2𝐇_\omega (𝐫).$$
(2)
In addition, if $`ϵ(𝐫)`$ and $`\mu (𝐫)`$ are periodic functions, following the derivation of Bloch theory, Eq. (2) can be changed to

$`{\displaystyle \underset{𝐆^{}}{}}\left(𝐤+𝐆\right)\times ϵ_{𝐆𝐆^{}}^1\left(𝐤+𝐆^{}\right)\times 𝐇_{𝐤,𝐆^{}}`$
$`=\omega ^2{\displaystyle \underset{𝐆^{}}{}}\mu _{𝐆𝐆^{}}𝐇_{𝐤,𝐆^{}},`$ (3)
where $`𝐇_\omega \left(𝐫\right)=\underset{𝐆}{}e^{i\left(𝐤+𝐆\right)𝐫}𝐇_{𝐤,𝐆}`$, $`ϵ\left(𝐫\right)=\underset{𝐆}{}e^{i𝐆𝐫}ϵ_𝐆`$, $`\mu \left(𝐫\right)=\underset{𝐆}{}e^{i\mathrm{𝐆𝐫}}\mu _𝐆`$, and $`\{𝐆\}`$ is the set of the reciprocal lattice. Since in this paper we consider only two-dimensional cases, we have $`k_z=0`$, and the electromagnetic waves can be decoupled as E polarization (TE) and H polarization (TM) modes. For example, the TM mode of $`𝐇`$ field is written as $`𝐇=H_z\widehat{z}`$ and satisfy

$`{\displaystyle \underset{𝐆^{}}{}}ϵ_{𝐆𝐆^{}}^1\left(𝐤+𝐆\right)\left(𝐤+𝐆^{}\right)H_{z,𝐤,𝐆^{}}`$
$`=\omega ^2{\displaystyle \underset{𝐆^{}}{}}\mu _{𝐆𝐆^{}}H_{z,𝐤,𝐆^{}}.`$ (4)
Conventionally, treating Eq. (4) as an eigenvalue problem, the propagating bulk modes of an infinite periodic system for a given real $`𝐤`$ can be obtained straightforwardly. However, in most of the situations we also need to know the transmittance and reflectance of a finite or semi-infinite system for a incident light. In order to obtain these quantities, various methods such like LKKR 8 ; 9 ; 11 , transfer matrix 12 ; 13 ; 14 ; 15 , and scattering matrix 16  method have been proposed. In these kinds of methods, a photonic crystal slab is treated as the stack of many gratings. The matrix problem for only one grating layer is solved first, then multiplying the matrices layer by layer and the total transmittance and reflectance can be determined. Although these methods are successful, however, in order to confirm the good numerical accuracy, the number of layers should be restricted to a small value. In addition, using these methods it is hard to find the relations between the transmittance and the original band structure.

In this paper, we propose an alternative method to calculate the transmitted and reflected waves from the interface between vacuum and a semi-infinite photonic crystal for a incident plane wave. The method is based on Eq. (4), which contains all the information of the band structure of the system. We thus can easily analyze the physical meanings of various phenomena using the knowledge getting from the solid state physics.

Here, if the $`k_x`$ and $`k_y`$ are fixed, the frequency $`\omega `$ can be solved as an eigenvalue problem. If the system is infinitely extended, there are just only propagation modes that can survive. However, sometimes we have to calculate the transmission and reflection coefficients for a finite sized or semi-infinite sample, which have at least one boundary. At the edges of the sample, the periodic structure is broken and the evanescent modes must be considered. However, it is impossible to obtain the evanescent solutions from the eigenvalue problem of Eq. (4) because it provides the solutions for an extended bulk and the boundary conditions at infinite restrict the $`k_x`$ and $`k_y`$ to be taken as real values.

Nevertheless, can we make some modifications to produce the attenuated solutions of Eq. (4)? The answer is yes. For different purposes, there are two kinds of calculations: one is to fix the direction of k and frequency; the other is to fix $`k_y`$ and frequency. With a simple transformation, the Eq. (4) can be rewritten as

$`\left(\begin{array}{cccccccccccccccccccc}0& \widehat{𝐈}& & & & & & & & & & & & & & & & & & \\ ϵ_{𝐆𝐆^{\prime \prime }}\left[\omega ^2\mu _{𝐆^{\prime \prime }𝐆^{}}ϵ_{𝐆^{\prime \prime }𝐆^{}}^1𝐆^{\prime \prime }𝐆^{}\right]& ϵ_{𝐆𝐆^{\prime \prime }}\left[ϵ_{𝐆^{\prime \prime }𝐆^{}}^1\widehat{𝐤}\left(𝐆^{\prime \prime }+𝐆^{}\right)\right]& & & & & & & & & & & & & & & & & & \end{array}\right)\left(\begin{array}{c}H_𝐆^{}\\ kH_𝐆^{}\end{array}\right)=k\left(\begin{array}{c}H_𝐆^{}\\ kH_𝐆^{}\end{array}\right)`$ (13)
for fixed direction and frequency, where (respectively) $`H_𝐆^{}`$, $`\widehat{𝐈}`$ and $`\widehat{𝐤}`$ denote the abbreviation of $`H_{z,𝐤,𝐆^{}}`$,$`\delta _{𝐆,𝐆^{}}`$ and unit vector of $`𝐤`$, and

$`\left(\begin{array}{cccccccccccccccccccc}0& \widehat{𝐈}& & & & & & & & & & & & & & & & & & \\ ϵ_{𝐆𝐆^{\prime \prime }}\left[\omega ^2\mu _{𝐆^{\prime \prime }𝐆^{}}ϵ_{𝐆^{\prime \prime }𝐆^{}}^1\left(𝐆^{\prime \prime }+k_y\widehat{𝐲}\right)\left(𝐆^{}+k_y\widehat{𝐲}\right)\right]& \widehat{𝐏}& & & & & & & & & & & & & & & & & & \end{array}\right)\left(\begin{array}{c}H_𝐆^{}\\ k_xH_𝐆^{}\end{array}\right)=k_x\left(\begin{array}{c}H_𝐆^{}\\ k_xH_𝐆^{}\end{array}\right)`$ (22)
for fixed $`k_y`$ and frequency, where $`\widehat{𝐲}`$ and $`\widehat{𝐏}`$ denote unit vector of y direction and $`ϵ_{𝐆𝐆^{\prime \prime }}\left[ϵ_{𝐆^{\prime \prime }𝐆^{}}^1\left(𝐆_x^{\prime \prime }+𝐆_x^{}\right)\right]`$, respectively. If the $`𝐤`$ and $`𝐆`$ in Eq. (13) are acted by a rotation operator $`\widehat{\mathrm{\Theta }}`$ which rotates $`\widehat{𝐤}`$ to x-direction — i.e. $`\widehat{\mathrm{\Theta }}\widehat{𝐤}=\widehat{𝐱}`$ — and we define $`\stackrel{~}{𝐆}\widehat{\mathrm{\Theta }}\widehat{𝐆}`$, then Eq. (13) becomes Eq. (22) in which $`\stackrel{~}{𝐆}`$ is substituted for $`𝐆`$ and $`k_y=0`$. Following this cue, the Eq. (13) can be considered as a master equation for solving a problem where the incident light is always perpendicular to the interface. From which we can easily determine the penetration depth along direction $`\widehat{𝐤}`$. If our purpose is to use the band gap effect of the structure, then the result obtained from this calculation will tell us how to cut the sample to get the highest performance.

On the other hand, Eq. (22) can be used to deal with the problem for light incident with different angle $`\theta `$, which is the angle between the normal vector of the interface and $`𝐤`$ of the incident light. Here $`\mathrm{tan}\theta =k_y/\sqrt{\omega ^2/c^2k_y^2}`$. Figure 1 explains the details.

Since we can obtain all the eigenvectors of the system, the transmission and reflection spectra can also be obtained. Based on the continuity conditions at the interface, the relationship between the $`𝐇`$ fields in region I and region II can be written as

$`\left(\begin{array}{cccccccccccccccccccc}x_0y|H_m^\text{I}& x_0y|H_m^{\text{II}}& & & & & & & & & & & & & & & & & & \\ x_0y|ϵ^1_x|H_m^\text{I}& x_0y|ϵ^1_x|H_m^{\text{II}}& & & & & & & & & & & & & & & & & & \end{array}\right)\left(\begin{array}{c}H_m^\text{I}|\widehat{𝐑}|H_0^\text{I}\\ H_m^{\text{II}}|\widehat{𝐓}|H_0^\text{I}\end{array}\right)=\left(\begin{array}{c}x_0y|H_0^\text{I}\\ x_0y|ϵ^1_x|H_0^\text{I}\end{array}\right),`$ (31)
where $`\widehat{𝐑}`$ and $`\widehat{𝐓}`$ are the reflection and transmission operators, $`H_m^\text{I}`$ and $`H_m^{\text{II}}`$ are the reflection modes in region I and transmission modes in region II, the $`m`$ denotes the different modes, respectively. And $`H_0^\text{I}`$ is the incident light.

If Eq. (31) is expanded in K-space, it can be rewritten as

$`\widehat{𝐀}\left(\begin{array}{c}H_{\text{I}m}|\widehat{𝐑}|H_{\mathrm{I}0}\\ H_{\text{II}m}|\widehat{𝐓}|H_{\mathrm{I}0}\end{array}\right)=\left(\begin{array}{c}H_{z,𝐤,G_y}^{\text{I},i}\\ k_{G_y,x}^{\text{I,}i}H_{z,𝐤,G_y}^{\text{I},i}\end{array}\right),`$ (38)
where

$$\widehat{𝐀}=\left(\begin{array}{cccccccccccccccccccc}H_{z,𝐤,G_y}^{\text{I},r}& \underset{G_x}{}H_{z,𝐤,𝐆}^{\text{II},t}& & & & & & & & & & & & & & & & & & \\ k_{mx}^\text{I}H_{z,𝐤,G_y}^{\text{I},r}& \underset{G_x,𝐆^{}}{}ϵ_{𝐆𝐆^{}}^1\left(k_{mx}^{\text{II}}+G_x^{}\right)H_{z,𝐤,𝐆}^{\text{II},t}& & & & & & & & & & & & & & & & & & \end{array}\right),$$

and $`H_{z,𝐤,𝐆}`$ can be gotten from Eq. (22), and $`i,r`$ and $`t`$ denote incident, reflection and transmission, respectively. To determine the transmission and reflection coefficients, we have to first decide the direction of the Poynting vector of every mode. For $`𝐤`$ is a real vector, it can either be obtained from $`𝐯_g=\stackrel{}{}_𝐤\omega `$ or from

$`{\displaystyle \underset{cell}{}}\mathrm{Re}\left\{{\displaystyle \frac{i}{\omega ϵ}}H_z^{}\stackrel{}{}H_z\right\}𝑑r^2`$
$`={\displaystyle \underset{𝐆,𝐆^{}}{}}\mathrm{Re}\left\{{\displaystyle \frac{1}{\omega }}H_{z,𝐤,𝐆}^{}ϵ_{𝐆𝐆^{}}^1\left(𝐤+𝐆^{}\right)H_{z,𝐤,𝐆^{}}\right\}`$ (39)
For $`𝐤`$ is a complex number, propagation toward right hand side corresponds to $`\mathrm{Im}(k_x)>0`$. When group velocity and $`H_{z,𝐤,𝐆}`$ of each mode are known, the transmittance $`𝐓`$ and reflectance $`𝐑`$ can be gotten from them, and the accuracy can be estimated from how $`𝐑+𝐓`$ is close to one. Sometimes, $`det|𝐀|`$ is possible to become zero, if so, the Eq. (38) has nonezero solutions when there is no incident light. This kind of wave is a surface state which resembles the surface plasmon propagating along the surface of a metal. However, we do not discuss it in this paper. One more interesting thing among these three equations — from Eq.(4) to Eq.(22) — is that they are identical to each other. Because the second row of matrices of left side of Eq. (13) and Eq. (22) are equal to Eq. (4) while the $`\{𝐆\}`$ in these equations are equal. Thus any eigenfunction of one of these three equations satisfies another two equations. It leads to two useful things: (i) the real k contours of these three methods with equal frequency are the same. (ii) The $`\{𝐆\}`$ needn’t to be changed when these three equations are considered as a series of policy tools. For example, we should decide where the band edge is and select the $`\omega _0`$ near the band edge to obtain the band structure from Eq. (4) if we hope to observe what happened near the band gap. By replacing the frequencies of Eq. (13) and (22) by $`\omega _0`$, the penetration depth can be derived from Eq. (13) and so do the variations of transmittance while the solutions of Eq. (22) are used in Eq. (31). During this process, even if we just select $`𝐆=0`$ and drop out $`𝐆0`$ for average assumption in the calculations of Eq. (4), the $`\{𝐆\}`$ still needn’t to be changed in the following calculations.

For simplicity, the structure we used in obtaining Fig. 2 and Fig. 3 is the square lattice with GaAs ($`ϵ=11.43`$ ) cylinders each with radius=$`0.15a`$ inside vacuum; whereas in Fig. 4 we employ vacuum cylinders each with radius $`0.43077a`$ inside the PbO($`ϵ=2.72`$ ) background. In all cases, $`a`$ is the lattice constant and the primitive vectors are given by $`𝐚_1=(a,0)`$ and $`𝐚_2=(0,a)`$.

The first application is a selection rule to determine the interface direction for the highest performance of light insulation. We solve Eq. (13) at a frequency equal to $`0.4(2\pi c/a)`$ and the results are presented in Fig. 3. In Fig. 2, there is a band gap at $`\omega =0.4(2\pi c/a)`$ for the TE mode, so that in Fig. 3(c), there are no real number solutions inside the first Brillouin zone. However, outside and far away from the 1st Brillouin zone such solutions exist, which are fake and are caused by the finite basis used in the calculations. To find the interface direction we first choose a direction $`\widehat{𝐤}`$ and use Eq. (13) to find a $`k`$ that has the smallest $`|\mathrm{Im}(k)|`$ value, which denoted as $`k_I(\theta )`$ and determines the main decay trend for a wave propagating along $`\widehat{𝐤}`$. The second step is to scan angles from $`0^{}`$ to $`45^{}`$ to find an angle $`\theta _0`$ that has the maximum $`k_I(\theta )`$. The details are shown in Fig. 3(d), where we calculate the TE modes, and the penetration depths (i.e., $`2\pi /k_I(\theta )`$) for $`0^{}`$ and $`45^{}`$ are $`6.8540a`$ and $`3.3272a`$, respectively. This indicates that when we produce a sample that cut along $`45^{}`$, it just needs 4 or 5 layers to stop the light with $`\omega =0.4(2\pi c/a)`$ for TE modes instead of 7 or 8 layers for $`0^{}`$.

The second application is to fix $`k_y`$ and frequency in order to obtain $`k_x`$. For comparison, we select the system discussed in 7  to contrast with our system and show the result in Fig. 4. The system in 7  is a 16 layers photonic crystal which is a square lattice (with a lattice constant $`a`$) of air columns (radius equals $`0.43077a`$ and is located at the center of a unit cell) in a dielectric substrate PbO placed in air and our system is a semi-infinite photonic crystal placed in air. By using Eq. (22) and Eq. (38) the transmittance can be obtained as shown in Fig. 4. Because our system is infinite, we can find something quite different.

(1) The solid lines and dashed lines are almost smooth curves except in the gap regime and at some special points ($`\omega =0.74`$ and $`0.85`$ in TE mode). The oscillating solid lines with dots in Fig. 4 represent the solutions of 7 . It is obvious that our curves are different from that of 7 . The oscillation behavior is owing to finite thickness. They can be easily explained by a roughly consideration of the average dielectric $`\overline{ϵ}ϵ_{cell}`$. For low frequency, the most important contribution of $`ϵ_𝐆`$ is $`ϵ_{𝐆=0}`$, which is equal to $`\overline{ϵ}`$, therefore the $`16`$ layers photonic crystals can be considered as an effective material with uniform dielectric $`\overline{ϵ}=1.7173ϵ_0`$ and width $`16a`$, where $`a`$ is the lattice constant of the photonic crystal and $`ϵ_0`$ is the permittivity of free space. This is a typical 1D problem and the waves in both sides of this material can be connected by transfer matrix whose formula is

$$𝐌_\omega \left(\begin{array}{cc}e^{ik16a}& 0\\ 0& e^{ik16a}\end{array}\right)𝐌_\omega ^1,$$

where $`𝐌_\omega `$ is a $`2\times 2`$ matrix with $`det|𝐌_\omega |0`$ determined by the optical impedance contrast and the incident angle, and $`k=\omega \sqrt{\overline{ϵ}\mu _0}`$ is the effective wave number, where $`\mu _0`$ is permeability of free space. It is obvious that the transfer matrices equals $`\pm 𝐈`$ when $`k16a=n\pi `$, where $`n`$ is an integer number. That means the frequency difference between two neighboring peaks is

$$\mathrm{\Delta }\stackrel{~}{\omega }=\frac{1}{32\sqrt{\overline{ϵ}/ϵ_0}}=\frac{1}{32\sqrt{1.7173}}0.0238,$$

where $`\stackrel{~}{\omega }=\frac{\omega a}{2\pi c}`$. Comparing with the average width $`0.0222`$ of peaks between $`\frac{\omega a}{2\pi c}=0.5`$ to $`0.7`$, we can say the oscillation in Fig. 4 comes from the effect of finite size.

(2) In the vicinity of $`\stackrel{~}{\omega }0.74`$, according to the band structure calculation results the spectrum should ascend rapidly when $`\stackrel{~}{\omega }`$ is increasing, because it is at band edge. But, the real situation appears in Fig. 4(a) is ascending quickly and descending immediately to near zero transmittance. Our explanation is that there are two propagation mode’s with group velocity $`v_g0`$, so they do not contribute to the transmittance.

(3) In Fig. 4(a), the valley of transmittance near $`\omega =0.85`$ disappears when the interface is chosen to pass through the center of the vacuum cylinders. This reveals that it is possible to stop the light at some isolated frequency points by appropriately choosing the cutting plane of the photonic crystal even if the frequencies are outside of band gaps.

(4) In Fig. 4(b), there should be a forbidden band for $`\omega =0.75`$ to 0.78, but the line with dots does not show this result. According our calculation the wave attenuation rate for a 16-layer structure is about 0.0733, which agrees with the result shown by the line with dots. This phenomenon shows that the evanescent modes do contribute the transmittance in a finite thickness structure.

In these two applications, the evanescent modes are necessary and useful for the calculations of semi-infinite system, and this method also provides us some information such like for how large a separation distance between two defects can they be treated as independent in super-cell method.

In conclusion, the method we present here may not be efficient enough, because in the calculation we get results both inside and outside of the first Brillouin zone (FBZ). However, only the results inside the FBZ are useful and the others repeat the same information and are redundant. For example, if we use $`N^2`$ bases, there are only $`2N`$ useful eigenfunctions. Besides, this method has several advantages. First, from this method we can easily realize and analyze some properties of periodic systems with interface and the computational time is independent of the number of layers. Thus, even if the number of layers is very large, it will save much time. Second, using Eq. (39), we can also calculate the density of states, D($`\omega `$), through

$`D\left(\omega \right)={\displaystyle \underset{shell}{}}{\displaystyle \frac{dk_{//}}{\left|_𝐤\omega \right|}}.`$ (40)
They are especially useful when we aim to calculate the density of states in some small frequency regimes.

We are now investigating the cases of finite size specimens and a structure with line defects.

Finally, we thank Prof. B. Y. Gu for instructing us about Andreev reflection, which gave us a chance to employ this idea about Eq. (13), and thank Dr. P. G. Luan who let us find more possibilities with this method.