text:
# I Introduction

## I Introduction

A number of pattern formation problems involving an interface are computationally difficult despite the fact that the bulk equations are linear. Examples include Laplacian dynamics (diffusion-limited aggregation, limit of vanishing Péclet number in dendritic growth), flow problems in the Stokes approximation, and elastic problems. They are difficult because of the long-range self-interaction of the interface mediated by the external field. With the availability of modern computer power, these problems can usually be solved to good accuracy in two dimensions. But none of them has been successfully treated on a large scale in three dimensions so far.

Whenever an interface equation can be derived, this offers the advantage of better tractability, both numerically and analytically. Interface equations for three-dimensional systems are two-dimensional, rendering much larger systems accessible than bulk simulations kassner98 . For problems containing a length scale restricting the interaction distance, such as the diffusion length in solidification problems, interface equations are local in the long-wave limit misbah91 ; kassner94 , i.e., they are partial differential equations. The Laplace or Lamé equations, however, do not have a cutoff distance. Hence, it is interesting and important to investigate how this feature affects the long-wave limit, if any. There have been interface equations containing nonlocal but linear terms in the literature before sivashinsky77 ; hobbs92 ; marietti01a . As it will turn out, the nonlocal aspects of the system we consider are much more ferocious, leading to terms that are both nonlocal and nonlinear. We are also aware of the case of a previously derived amplitude equation containing nonlocal nonlinearities elmer88 ; christen90 . However, the equation in question is weakly nonlinear whereas we will obtain a strongly nonlinear equation.

Let us recall that close to the instability threshold for the merging of order, systems exhibiting a type Is instability in the nomenclature of Cross and Hohenberg cross93 , i.e., an instability at finite wave number, are all expected to be described by a universal amplitude equation of the Ginzburg-Landau type. Otherwise, when the critical wave number approaches zero (type IIs) the dynamics is described by a long-wave equation such as the Kuramoto-Sivashinsky (KS) equation. A prerequisite for these derivations is that nonlocal interactions are screened off beyond some length (small in comparison with $`1/q_c`$). When long-range interactions (e.g. in electrostatic, magnetic, elastic systems, etc.) are effectively present, no universal equation has been known as yet for the latter case. The results of elmer88 ; christen90  pertain to the former case, i.e., they constitute a generalization of the Ginzburg-Landau description including nonlocal interactions. We derive here a universal form of dynamics in the long-wavelength limit when long-range interactions prevail and obtain a result that is complementary to the Kuramoto-Sivashinsky description. The derivation is exemplified on a system undergoing a surface instability due to elastic stress. Besides the universal feature of the resolved question, we present the practical virtue of this strategy.

Applying our approach to a particular system, a uniaxially strained solid undergoing the Grinfeld instability grinfeld86 ; torii92 , a problem which has recently attracted the community of crystal growth nozieres93 ; mueller99 ; kassner99 , we obtain an equation that, when truncated to the leading-order nonlinearities should provide an example of the universal dynamics. As we shall see, this equation contains fewer parameters than to be expected which is most likely due to the fact that in the case of the Grinfeld instability the dynamics is variational. Introducing an additional free parameter, we are able to give the generic form of the universal equation. A derivation of this equation from symmetry considerations will be presented elsewhere kassner02 .

The paper is organized as follows. In Sec. II, we write down the model equations for the elastically strained system and introduce the appropriate rescaling to perform an asymptotic analysis. Section III describes the asymptotic expansion and matching procedure and presents the final interface equation for two dynamical situations corresponding to conservative and to nonconservative dynamics. In Sec. IV, we give some simulation results for these equations. Section V contains a discussion of universality aspects and gives an example of a simulation of the generic equation. Finally, we will briefly summarize some conclusions and give an outlook in Sec. VI.

## II  Model equations and nondimensionalization

The physics of the Grinfeld instability has been discussed in some detail in durand98 . We consider the case of a solid in contact with its melt; let us call the interface position $`\zeta (x,t)`$ (i.e., we restrict ourselves to two dimensions). It is related to the stress distribution via nozieres93 ; durand96 ; durand98 ; kassner01 :

$`{\displaystyle \frac{\zeta _t(x,t)}{\sqrt{1+\zeta _x^2}}}`$ $`=`$ $`{\displaystyle \frac{1}{k\rho _s}}\left({\displaystyle \frac{1\nu ^2}{2E}}\left[(\sigma _{tt}\sigma _{nn})^2\sigma _0^2\right]+\gamma \kappa +\mathrm{\Delta }\rho g\zeta \right).`$ (1)
Herein, $`1/k`$ is a mobility, $`\rho _s`$ the density of the solid, $`E`$ is Young’s modulus, $`\nu `$ the Poisson number, $`\gamma `$ is the surface tension, $`\mathrm{\Delta }\rho `$ the density contrast between solid and liquid, $`g`$ the gravitational acceleration. $`\kappa =\zeta _{xx}/\left(1+\zeta _x^2\right)^{3/2}`$ is the interface curvature, and partial derivatives of $`\zeta `$ are denoted by subscripts. $`\sigma _{tt}`$ and $`\sigma _{nn}`$ are tangential and normal stresses at the interface, i.e., $`\sigma _{nn}=n_i\sigma _{ij}n_j`$, $`\sigma _{tt}=t_i\sigma _{ij}t_j`$, where $`n_i`$ and $`t_i`$ are the components of the normal and the tangent to the interface. By $`\sigma _0`$, we denote the externally imposed uniaxial stress, and the interface position is measured from the equilibrium position of a planar interface. (Setting $`\kappa =0`$ and $`\left|\sigma _{tt}\sigma _{nn}\right|=\left|\sigma _0\right|`$ in (1) leads to a steady state $`\zeta =0`$.) Equation (1) describes the mass current density which is proportional to chemical potential changes produced by strain, capillarity and gravity.

To solve this equation, it must be supplemented with the equations describing mechanical equilibrium. Assuming linear isotropic elasticity,

$$\sigma _{ij}=\frac{E}{1+\nu }\left(u_{ij}+\frac{\nu }{12\nu }u_{kk}\delta _{ij}\right),$$
(2)
where $`u_{ij}`$ is the strain tensor, expressible via the elastic displacements $`u_i`$ according to $`u_{ij}=\frac{1}{2}\left(_iu_j+_ju_i\right),`$ the condition for mechanical equilibrium, $`_j\sigma _{ij}=0`$ (neglecting gravity in the bulk as a small cross effect), takes the form of the Lamé equations:

$$(12\nu )^2𝐮+(𝐮)=0.$$
(3)
Note that the static version is sufficient since all interface motions are slow in comparison with sound propagation. Boundary conditions are

$$\sigma _{nn}=\sigma _{nt}=0$$
(4)
at the interface, periodicity in the $`x`$ direction with some prescribed wavelength $`2\pi /q`$, and

$$\sigma _{ij}\sigma _0\delta _{ix}\delta _{jx}\text{for}z\mathrm{}.$$
(5)
Herein, $`\sigma _{nt}=n_i\sigma _{ij}t_j`$, is the shear stress at the interface, which must evidently vanish if the latter is in contact with a fluid. Another cross effect has been neglected here, viz. that of capillary overpressure, leading to a curvature dependent jump of $`\sigma _{nn}`$ across the interface. An overall constant pressure, not influencing the dynamics of the interface has been subtracted out of the definition of $`\sigma _{nn}`$.

Due to the free boundary character, equations (2) through (5) do not completely specify the problem and there is a need for an additional condition, which is just equation (1) expressing the fact that the chemical potential difference between solid and liquid is proportional to the normal mass current across the interface. This means that we assume the interface to be microscopically rough (so that a linear relation holds between the current and its conjugate variable, the chemical potential difference) and the solid to be in contact with its liquid or vapor.

Note however, that for the solution of the elastic problem at one instant of time with given interface position we do not need (1) (since we have neglected sound propagation effects). This suggests to attack the problem in two stages: first solve the elastic equations, then move the interface according to the elastic fields, and repeat the procedure for the new interface position.

The question is how to obtain a closed interface equation, and to identify its generic form. In order to proceed, we first nondimensionalize the above equations by introducing two length scales $`\mathrm{}_1=\gamma E/\sigma _0^2(1\nu ^2)`$ (the Griffith length griffith20 , apart from a factor of $`4/\pi `$) and $`\mathrm{}_2=(1\nu ^2)\sigma _0^2/2\mathrm{\Delta }\rho gE`$ (a gravity length) as well as a time scale $`\tau =k\rho _s\mathrm{}_1^2/\gamma `$ kassner01 . Then we use Hooke’s law (2) to express stresses by strains and reduce the latter by a dimensionless factor $`E/(1+\nu )\sigma _0`$. We find

$$\frac{\zeta _t(x,t)}{\sqrt{1+\zeta _x^2}}=\left\{\frac{1}{2}\left[\left(u_{tt}u_{nn}\right)^21\right]+\kappa +\alpha \zeta \right\},$$
(6)
where $`\alpha =\mathrm{}_1/2\mathrm{}_2`$, i.e. $`\alpha 0`$. A linear stability analysis of the full set of equations (2) through (6) yields the dispersion relation kassner01 :

$$\omega =2|q|q^2\alpha ,$$
(7)
where $`\omega `$ and $`q`$ are the nondimensional growth rate and wavenumber, respectively. The absolute value of $`q`$ arises, because the eigensolutions of the elastic problem carry a factor $`\mathrm{exp}(\pm iqx+|q|z)`$, approaching zero for $`z\mathrm{}`$. A linear instability will arise for $`\alpha <1`$.

As is obvious from the derivation, to be able to manifest itself the instability requires mass transport. Besides the case described of a liquid in contact with its melt, corresponding to nonconserved dynamics, another situation is of particular interest, due to possible applications in epitaxial growth. This is a solid in vacuum. Material transport will ordinarily be dominated by surface diffusion then. Consequently, we have conserved dynamics, with eq. (6) replaced by

$$\frac{\zeta _t(x,t)}{\sqrt{1+\zeta _x^2}}=\frac{1}{\sqrt{1+\zeta _x^2}}_x\frac{1}{\sqrt{1+\zeta _x^2}}_x\left\{\frac{1}{2}\left[\left(u_{tt}u_{nn}\right)^21\right]+\kappa +\alpha \zeta \right\},$$
(8)
and the dispersion relation (7) is simply multiplied by a factor $`q^2`$ on the right-hand side. Since the gravity term is normally negligible in experiments on epitaxial growth (unless a temperature gradient mimicks a strong gravitational field durand96 ; durand98 ), we may set $`\alpha =0`$ in this case, obtaining a parameter-free equation.

As it turns out, the asymptotic analysis can be performed together for both dynamical situations, since its most extensive part consists in solving the elastic problem (which is the same for both cases). Only afterwards is the solution inserted into the equation describing the interface evolution. At intermediate steps of the calculation, we will therefore write down equations for the nonconservative case only, but we shall give the final interface equation for both cases.

Let us now consider a situation, where $`\mathrm{}_1`$ is very large, i.e., we introduce a small parameter $`ϵ`$, setting $`\mathrm{}_1=O(1/ϵ)`$. This corresponds to a very small stress, but we can nevertheless reach an unstable state by reducing gravity or the density difference, so $`\mathrm{}_2`$ becomes larger than $`\mathrm{}_1/2`$.

Then it will be natural to measure length scales in units of $`ϵ\mathrm{}_1`$ \[which is $`O(1)`$\] rather than $`\mathrm{}_1`$ and time scales in units of $`ϵ^2\tau `$. Coordinates are transformed according to $`\stackrel{~}{x}=x/ϵ`$, $`\stackrel{~}{z}=z/ϵ`$, and $`\stackrel{~}{t}=t/ϵ^2`$. Referring to the new coordinates, the dispersion relation transforms into $`\stackrel{~}{\omega }=2ϵ|\stackrel{~}{q}|\stackrel{~}{q}^2ϵ^2\alpha `$. For small $`ϵ`$, each term in this relation must behave as $`ϵ^2`$, hence we have for the wavenumber $`\stackrel{~}{q}ϵ`$. A long-wave equation should therefore be derivable. To perform the asymptotics, we introduce another set of coordinates via $`X=ϵ\stackrel{~}{x}`$, $`Z=\stackrel{~}{z}\zeta (x,t)`$, $`T=ϵ^2\stackrel{~}{t}`$, conveniently mapping the interface position to $`Z=0`$. $`X`$ and $`T`$ are slow variables. We set $`U=u_x(1\nu )X`$, $`V=u_z+\nu (Z+\zeta )`$, where $`u_x`$ and $`u_z`$ are the displacements, which makes the boundary conditions at $`\mathrm{}`$ homogeneous ($`U0`$, $`V0`$). The transformed equations (1,3) read:

$`ϵ\zeta _T(1+ϵ^2\zeta _X^2)^{3/2}=`$
$`\{{\displaystyle \frac{1}{2}}[(1ϵ^2\zeta _X^2)(1+_XU\zeta _X_ZU_ZV)+2\zeta _X(_ZU+ϵ^2(_XV\zeta _X_ZV))]^2`$
$`{\displaystyle \frac{1}{2}}(1+ϵ^2\zeta _X^2)^2ϵ\zeta _{XX}(1+ϵ^2\zeta _X^2)^{1/2}+ϵ\alpha \zeta (1+ϵ^2\zeta _X^2)^2\},`$ (9)
$`(12\nu )_Z^2U+ϵ^2\left[2(1\nu )\left(_X\zeta _X_Z\right)^2U+\left(_X\zeta _X_Z\right)_ZV\right]`$ $`=`$ $`0,`$ (10)
$`2(1\nu )_Z^2V+\left(_X\zeta _X_Z\right)_ZU+ϵ^2(12\nu )\left(_X\zeta _X_Z\right)^2V`$ $`=`$ $`0,`$ (11)
and the boundary conditions at the interface ($`Z=0`$) become

$`0`$ $`=`$ $`(1\nu )_ZV+\nu _XU(1\nu )\zeta _X_ZU`$ (12)
$`+ϵ^2\left[(12\nu )\left(\zeta _X^2\zeta _X_XV\right)+(1\nu )\zeta _X^2\left(_ZV+_XU\zeta _X_ZU\right)\right],`$
$`0`$ $`=`$ $`\left[_ZU+ϵ^2\left(_XV\zeta _X_ZV\right)\right]\left(1ϵ^2\zeta _X^2\right)+2ϵ^2\zeta _X\left(_ZV_XU+\zeta _X_ZU1\right).`$ (13)
Equations (9) through (13) constitute our starting point for the asymptotic analysis to be carried out.

## III  Asymptotic analysis

These equations are singular in $`Z`$, i.e., the boundary conditions at infinity, where the implicit assumption $`Z=O(1)`$ does not hold anymore, can be satisfied only with the trivial solution $`\zeta _X=0`$. This difficulty usually does not appear with equations that have an internal length scale, in contrast to the elastic equations which are devoid of such a scale. To overcome it, one introduces an inner region, where $`Z=O(1)`$, and an outer region, where $`\eta ϵZ=O(1)`$ hobbs92 , attempting then to match the solutions in their common domain of validity. In the inner region, Eqs. (10)–(13) hold, and solving the bulk equations may be reduced to ordinary differential equations in $`Z`$ via expansion in powers of $`ϵ`$. In the outer region, derivatives with respect to $`\eta `$ and $`X`$ are the same order of magnitude. Therefore, true partial differential equations have to be solved. It is here that the linearity of the Lamé equations becomes important, because it allows one to solve these equations in terms of Fourier transforms. Nontrivial aspects of the calculation then are largely due to the asymptotic matching of real space functions with Fourier transforms.

Denoting field variables in the outer region by small letters, $`u(X,\eta ,T)=U(X,Z,T)`$, $`v(X,\eta ,T)=V(X,Z,T)`$, we obtain the outer equations

$`(12\nu )_\eta ^2u+2(1\nu )\left(_Xϵ\zeta _X_\eta \right)^2u+ϵ\left(_Xϵ\zeta _X_\eta \right)_\eta v`$ $`=`$ $`0,`$ (14)
$`\left(_Xϵ\zeta _X_\eta \right)_\eta u+2(1\nu )ϵ_\eta ^2v+ϵ(12\nu )\left(_Xϵ\zeta _X_\eta \right)^2v`$ $`=`$ $`0,`$ (15)
the inner ones being given by (10,11). Expanding the interface position and the fields in powers of $`ϵ`$, we find that the expansion is singular in another sense. If it is assumed that $`\zeta =O(1)`$, nonlinear terms will not arise in the interface equation, a fact that can also be inferred from symmetry considerations and power counting. Scalings that lead to a nontrivial result are $`\zeta =O(ϵ^1)`$, $`u=O(1)`$, $`v=O(ϵ^1)`$. We therefore write

$$\zeta (X,T)=ϵ^1\zeta _1(X,T)+\zeta _0(X,T)+ϵ\zeta _1(X,T)+\mathrm{}$$
(16)
implying $`_{\stackrel{~}{x}}\zeta (X,T)=O(1)`$. Next we set

$`V(Z)`$ $`=`$ $`ϵ^1V_1(Z)+V_0(Z)+ϵV_1(Z)+\mathrm{}=v(\eta )=ϵ^1v_1(\eta )+v_0(\eta )+ϵv_1(\eta )+\mathrm{},`$
$`U(Z)`$ $`=`$ $`U_0(Z)+ϵU_1(Z)+\mathrm{}=u(\eta )=u_0(\eta )+ϵu_1(\eta )+\mathrm{},`$ (18)
where for brevity we have suppressed the $`X`$ and $`T`$ dependences. We then obtain as matching conditions (taking into account the $`ϵ`$ dependence of $`\eta `$)

$`U_0(Z)`$ $``$ $`u_0(0),(Z\mathrm{}),`$ (19)
$`V_1(Z)`$ $``$ $`v_1(0),(Z\mathrm{}),`$ (20)
$`U_1(Z)`$ $``$ $`Zu_0^{}(0)+u_1(0),(Z\mathrm{}),`$ (21)
$`V_0(Z)`$ $``$ $`Zv_1^{}(0)+v_0(0),(Z\mathrm{}),`$ (22)
where the prime denotes a derivative with respect to the argument ($`\eta `$). To the two leading orders, the inner equations are solved by

$`U_0(X,Z,T)`$ $`=`$ $`B_0(X,T),`$ (23)
$`V_1(X,Z,T)`$ $`=`$ $`D_1(X,T),`$ (24)
$`U_1(X,Z,T)`$ $`=`$ $`ZA_1(X,T)+B_1(X,T),`$ (25)
$`V_0(X,Z,T)`$ $`=`$ $`ZC_0(X,T)+D_0(X,T).`$ (26)
The outer equations need to be solved only at lowest order in $`ϵ`$. For brevity, we rename $`\zeta _1(X,T)`$ back to $`\zeta (X,T)`$. Confusion should not arise, since we will never need $`\zeta _0`$ nor $`\zeta _1`$. It should be kept in mind, however, that the new $`\zeta (X,T)`$ is $`O(1)`$, because it is the prefactor of $`ϵ^1`$. The outer solution then reads

$`u_0(X,\eta ,T)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[a_0(q,T)+b_0(q,T)(\eta +\zeta (X,T))\right]e^{\mathrm{i}qX+\left|q\right|(\eta +\zeta (X,T))}𝑑q,`$ (27)
$`v_1(X,\eta ,T)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left\{\mathrm{i}\mathrm{sign}(q)\left[a_0(q,T)+b_0(q,T)(\eta +\zeta (X,T))\right]+\mathrm{i}{\displaystyle \frac{34\nu }{q}}b_0(q,T)\right\}`$ (28)
$`e^{\mathrm{i}qX+\left|q\right|(\eta +\zeta (X,T))}dq.`$
Using the boundary conditions (12,13), we can express $`A_1`$ and $`C_0`$ by $`B_0`$ and $`D_1`$:

$`A_1(X)={\displaystyle \frac{1}{(1\nu )\left(1+\zeta _X^2\right)^2}}`$ $`\{\zeta _X[2(1\nu )+(12\nu )\zeta _X^2+((1\nu )\zeta _X^2+2\nu )B_0^{}(X)]`$ (29)
$`(1\nu \nu \zeta _X^2)D_1^{}(X)\},`$
$`C_0(X)={\displaystyle \frac{1}{(1\nu )\left(1+\zeta _X^2\right)^2}}`$ $`\left\{\zeta _X^2+\left(\nu +(1\nu )\zeta _X^2\right)\left(B_0^{}(X)+\zeta _X(X)D_1^{}(X)\right)\right\},`$ (30)
where we have left out the $`T`$ dependence (which is “passively” present everywhere but important only in the final equation – the Lamé equations are always solved at a fixed time).

The matching conditions at leading order give four equations expressing $`B_0`$, $`D_1`$, $`A_1`$, and $`C_0`$ by Fourier integrals involving $`a_0(q,T)`$ and $`b_0(q,T)`$. For $`B_0`$ and $`D_1`$, these integrals are essentially shown in Eqs. (27,28), we just need to set $`\eta =0`$ there. The two remaining equations read:

$`A_1(X)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[a_0(q)\left|q\right|+b_0(q)(1+\left|q\right|\zeta (X))\right]e^{\mathrm{i}qX+\left|q\right|\zeta (X)}𝑑q,`$ (31)
$`C_0(X)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\left\{\mathrm{i}qa_0(q)\mathrm{i}b_0(q)[\mathrm{sign}(q)+q\zeta (X)]+\mathrm{i}\mathrm{sign}(\mathrm{q})(34\nu )b_0(q)\right\}`$ (32)
$`e^{\mathrm{i}qX+\left|q\right|\zeta (X)}dq.`$
So we have six equations altogether, viz. Eqs. (27) through (32), for the six unknown functions $`B_0`$, $`D_1`$, $`A_1`$, $`C_0`$, $`a_0`$, and $`b_0`$ remark01 . Their solution determines the elastic fields, which then can be inserted into the interface equation (9). However, what we need is an explicit analytic solution. Its details are technical, so we just give some essential steps. First, we rewrite the equations in terms of Laplace transforms (rather: one-sided Fourier transforms)

$`\widehat{a}_0(x)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}a_0(q)e^{\mathrm{i}qx}𝑑q,`$ (33)
$`\widehat{b}_0(x)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{\mathrm{}}}b_0(q)e^{\mathrm{i}qx}𝑑q,`$ (34)
for example,

$$B_0(X)=\widehat{a}_0(X\mathrm{i}\zeta (X))+\zeta (X)\widehat{b}_0(X\mathrm{i}\zeta (X))+\mathrm{c}.\mathrm{c}..$$
(35)
Setting $`w=X+\mathrm{i}Z`$,

$`2\widehat{a}_0^{}(w)`$ $``$ $`m(X,Z)+\mathrm{i}n(X,Z),`$ (36)
$`2\widehat{b}_0(w)`$ $``$ $`o(X,Z)+\mathrm{i}p(X,Z),`$ (37)
with $`m`$, $`n`$, $`o`$, $`p`$ real functions, we can then reduce our problem to just two (real) equations involving complex quantities. These are the Eqs. (29,30) with all quantities expressed by $`\widehat{a}_0^{}(w)`$ and by $`\widehat{b}_0(w)`$, i.e., finally by $`m`$, $`n`$, $`o`$, and $`p`$.

Introducing the abbreviations

$`\lambda _1(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _X\left(22\nu +\left(12\nu \right)\zeta _{X}^{}{}_{}{}^{2}\right)}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (38)
$`\lambda _2(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _X\left(2\nu +\left(1\nu \right)\zeta _{X}^{}{}_{}{}^{2}\right)}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (39)
$`\lambda _3(X)`$ $`=`$ $`{\displaystyle \frac{1\nu \nu \zeta _{X}^{}{}_{}{}^{2}}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (40)
$`\lambda _4(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _{X}^{}{}_{}{}^{2}}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (41)
$`\lambda _5(X)`$ $`=`$ $`{\displaystyle \frac{\nu +\left(1\nu \right)\zeta _{X}^{}{}_{}{}^{2}}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (42)
$`\lambda _6(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _X\left(\nu +\left(1\nu \right)\zeta _{X}^{}{}_{}{}^{2}\right)}{\left(1\nu \right)\left(1+\zeta _{X}^{}{}_{}{}^{2}\right)^2}},`$ (43)
these two equations become

$`\left[1\lambda _2(X)\zeta _X+(34\nu )\lambda _3(X)\right]o(X,\zeta )+\left[\lambda _2(X)+\lambda _3(X)\zeta _X\right]m(X,\zeta )`$
$`=\lambda _1(X)+\left(\lambda _2(X)\lambda _3(X)\zeta _X\right)\zeta o_X(X,\zeta )`$
$`\left(1+\lambda _2(X)\zeta _X+\lambda _3(X)\right)\left(n(X,\zeta )+\zeta p_X(X,\zeta )\right)2(12\nu )\lambda _3(X)\zeta _Xp(X,\zeta ),`$
(44)
$`\left[1\lambda _5(X)+\lambda _6(X)\zeta _X\right]m(X,\zeta )+\left[\lambda _5(X)\zeta _X+(34\nu )\lambda _6(X)\right]o(X,\zeta )`$
$`=\lambda _4(X)\left(1\lambda _5(X)+\lambda _6(X)\zeta _X\right)\zeta o_X(X,\zeta )+2(12\nu )\left(1\lambda _6(X)\zeta _X\right)p(X,\zeta )`$
$`+\left(\lambda _5(X)\zeta _X+\lambda _6(X)\right)\left(n(X,\zeta )+\zeta p_X(X,\zeta )\right),`$ (45)
where $`p_X`$ and $`o_X`$ denote the derivatives of $`p(X,\zeta )`$, $`o(X,\zeta )`$ with respect to $`X`$ at fixed $`\zeta `$. That is, we first take the derivative and only then allow $`\zeta =\zeta (X)`$.

Since the four functions involved are real and imaginary parts of of analytic functions, only two of them are independent. Indeed, we have

$`n(X,\zeta )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{(X^{}X)m(X^{},0)}{(X^{}X)^2+\zeta ^2}}𝑑X^{},`$ (46)
$`p(X,\zeta )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{(X^{}X)o(X^{},0)}{(X^{}X)^2+\zeta ^2}}𝑑X^{}.`$ (47)
Equations (44) and (45), being linear, can be solved formally for $`o(X,\zeta )`$ and $`m(X,\zeta )`$. This leads to

$`o(X,\zeta )`$ $`=`$ $`r_1(X)\left[n(X,\zeta )+\zeta p_X(X,\zeta )\right]+r_2(X)p(X,\zeta )+r_3(X),`$ (48)
$`m(X,\zeta )`$ $`=`$ $`r_4(X)\left[n(X,\zeta )+\zeta p_X(X,\zeta )\right]+r_5(X)p(X,\zeta )+r_6(X)\zeta o_X(X,\zeta )`$ (49)
$`+r_7(X),`$
where the functions $`r_1`$ through $`r_7`$ are given by

$`r_1(X)`$ $`=`$ $`{\displaystyle \frac{1}{12\nu }},`$ (50)
$`r_2(X)`$ $`=`$ $`{\displaystyle \frac{2\zeta _X}{(12\nu )(1+\zeta _X^2)}},`$ (51)
$`r_3(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _X}{(12\nu )(1+\zeta _X^2)}},`$ (52)
$`r_4(X)`$ $`=`$ $`0,`$ (53)
$`r_5(X)`$ $`=`$ $`2{\displaystyle \frac{1+\nu (1+\zeta _{X}^{}{}_{}{}^{2})}{1+\zeta _{X}^{}{}_{}{}^{2}}},`$ (54)
$`r_6(X)`$ $`=`$ $`1,`$ (55)
$`r_7(X)`$ $`=`$ $`{\displaystyle \frac{\zeta _{X}^{}{}_{}{}^{2}}{1+\zeta _{X}^{}{}_{}{}^{2}}}.`$ (56)
The computations leading from (27) through (32) to (48,49) are straightforward but heavy. Some intermediate steps (e.g., the simplification of the $`r_i(X)`$) as well as a verification of the whole calculation at the end have been done using computer algebra (MAPLE). Due to the chosen scalings (16) – (18), Eqs. (48,49) constitute an exact reformulation of the outer elastic problem, which is linear in $`u`$ and $`v`$, so higher order equations have the same form as the leading order in $`ϵ`$. To obtain analytic results, we have to employ some approximation. We assume that $`\zeta (X)`$, the expansion coefficient of $`ϵ^1`$, is itself a numerically small quantity. Then expressions such as (46,47) become Hilbert transforms. Denoting the Hilbert transform of $`f(X)`$ by

$$\left[f(X)\right]\frac{1}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{f(X^{})}{X^{}X}dX^{},$$
(57)
where the integral is to be taken as a principal value, as indicated by the bar, we find, as the limit $`\zeta 1`$ of (48,49)

$`o(X,0)`$ $`=`$ $`r_1(X)\left[m(X,0)\right]r_2(X)\left[o(X,0)\right]+r_3(X),`$ (58)
$`m(X,0)`$ $`=`$ $`r_5(X)\left[o(X,0)\right]+r_7(X).`$ (59)
This set of linear equations can be solved iteratively. Assuming, and this is the second approximation, that $`\zeta _X(X)`$ is also numerically small, we may truncate the interation after the first step, taking accurately into account only terms up to order $`\zeta _X^2(X)`$. Using $`\left[[f(X)]\right]=f(X)`$, we can solve the truncated iteration to obtain explicit expressions for $`o(X,0)`$, $`m(X,0)`$ that read:

$`o(X,0)`$ $`=`$ $`(12\nu )r_3(X)+\left[r_7(X)\right],`$ (60)
$`m(X,0)`$ $`=`$ $`(12\nu )r_7(X)2(1\nu )(12\nu )\left[r_3(X)\right].`$ (61)
Substituting $`o`$, $`m`$, $`n=[m]`$, and $`p=[o]`$ back into the equations for $`A_1`$, $`C_0`$, $`B_0`$, and $`D_1`$, inserting these into (23) – (26) and then back into (9), we finally obtain the sought-for interface equation, in which we rename $`X`$ into $`x`$ and $`T`$ into $`t`$ (they are the same variables anyway):

$`{\displaystyle \frac{\zeta _t(x,t)}{(1+\zeta _x^2)^{1/2}}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{1{\displaystyle \frac{1}{(1+\zeta _x^2)^2}}\left[(1\zeta _x^2)\left(1+{\displaystyle \frac{2}{\pi }}\underset{\mathrm{}}{\overset{\mathrm{}}{{\displaystyle }}}{\displaystyle \frac{\zeta _x(x^{})}{(x^{}x)(1+\zeta _x(x^{})^2)}}dx^{}\right)+{\displaystyle \frac{4\zeta _x^2}{1+\zeta _x^2}}\right]^2\right\}`$ (62)
$`+{\displaystyle \frac{\zeta _{xx}}{(1+\zeta _x^2)^{3/2}}}\alpha \zeta .`$
The arguments of $`\zeta `$ ($`x`$ or $`x^{}`$ and $`t`$) have been suppressed wherever possible without ambiguity.

The nonlocal term drives the instability. It is interesting to see whether (62) will reproduce a coarsening dynamics of the type discussed in kassner01 . For in that case whether a groove grows or shrinks does not depend on local quantities (such as the curvature) alone. Hence, the scenario described cannot be expected from a merely local equation.

This is the melt-crystallization case kassner01 . It is easily checked that linearization of this equation produces the linear dispersion relation (7). The corresponding equation with dynamics controlled by surface diffusion is obtained simply by preceding the right-hand side with the operator $`\left[(1+\zeta _x^2)^{1/2}_x\right]^2`$:

$`\zeta _t(x,t)`$ $`=`$ $`_x{\displaystyle \frac{1}{\sqrt{1+\zeta _x^2}}}_x({\displaystyle \frac{1}{2}}\{{\displaystyle \frac{1}{(1+\zeta _x^2)^2}}[(1\zeta _x^2)(1+{\displaystyle \frac{2}{\pi }}\underset{\mathrm{}}{\overset{\mathrm{}}{{\displaystyle }}}{\displaystyle \frac{\zeta _x(x^{})}{(x^{}x)(1+\zeta _x(x^{})^2)}}dx^{})`$ (63)
$`+{\displaystyle \frac{4\zeta _x^2}{1+\zeta _x^2}}]^2\}{\displaystyle \frac{\zeta _{xx}}{(1+\zeta _x^2)^{3/2}}}+\alpha \zeta ).`$
For the system described by this equation, coarsening dynamics does not yet seem to have been studied in the literature, therefore it will be briefly discussed in the next section.

## IV  Some simulation results

Both equations (62) and (63) have been simulated numerically. Figure 1 gives an example for period-doubling dynamics found in the nonconservative dynamics when an initially periodic sinusoidal structure is submitted to a small perturbation of twice its basic periodicity.

A similar result is obtained for conservative dynamics and the same parameter $`\alpha `$, see Fig. 2. In both cases, $`\alpha `$ was close to 1, hence the instability was weak. For smaller values of $`\alpha `$, we found stronger differences between the behavior of the two systems. The main statement that can be made, however, is that in both of them similar coarsening scenarios are operative, since coarsening allows a reduction of the overall elastic energy. Because of the restriction that the average interface height has to remain constant in the conservative case, the latter has coarsening proceed more slowly, as can be verified by comparison of Figs. 1 and 2.

In Fig. 3, we give the dynamics obtained from a local perturbation of a periodic interface. The process of coarsening is more complex than the simple period-doubling scenario suggested in kassner01 . Frustration effects come into play when grooves compete via the nonlocal interaction, so it may happen that not the next-nearest neighbor of the leading groove is a winner of the competition but rather the third neighbor.

An extensive overview of the properties of Eqs. (62,63) is not intended here. It will be the subject of a different publication kassner02b . Suffice it to say that we find these amplitude equations to reproduce several qualitative features of the dynamics of the system faithfully, including nonlocally induced coarsening. Differences with the full dynamics kassner01  appear at large amplitudes. For vanishing gravity (i.e., $`\alpha =0`$), the amplitude is not predicted to saturate in the full model but does so in the amplitude equation. This is not so surprising given that the latter was derived assuming small amplitudes (more precisely a small prefactor of the $`O(ϵ^1)`$ term in the amplitude expansion). In principle, systematic improvement is possible by continuing the iteration scheme truncated after the first step and including higher powers of $`\zeta `$ and $`\zeta _x`$. In practice, this improvement results in additional technical difficulties.

It is certainly not possible to claim genericity of Eqs. (62,63). Their derivation accounts for only the lowest-order (i.e., quadratic) nonlinearities exactly, while resumming some nonlinearities to infinite order, which are presumably specific for the elastic system. However, if there is a universal amplitude equation for the systems considered, then it must contain the terms that our derivation produces exactly, and it may be sufficient to consider just these terms.

## V  Generic equations

Before looking at the generic equation(s) derivable from our particular case, let us briefly recapitulate the general situation concerning amplitude equations of universal validity.

There is an over-abundance of nonequilibrium systems which spontaneously build up an organized pattern from an initially structureless state when they are taken sufficiently far from equilibrium cross93 . Typical examples are present in hydrodynamical systems (e.g., Rayleigh-Bénard convection), chemical reactions (e.g., Turing systems), crystal growth (the Mullins-Sekerka instability), and so on. Despite the fact that the underlying physical and chemical mechanisms are diverse in these systems, sufficiently close to the instability threshold they all fall within the same universality class. Let $`q_c`$ denote the wavenumber of the emerging ordered pattern, and $`h(x,t)`$ a field describing the pattern (which can stand for a component of the velocity field in hydrodynamics, or an interface position in crystal growth, etc…) along the $`x`$ direction. Close to the threshold it is known cross93  that $`hA(x,t)e^{iq_cx}`$, where $`e^{iq_cx}`$ describes the (rapid) periodic variation of the field due to the birth of order, while $`A`$ is a slowly varying (slow with respect to $`e^{iq_cx}`$) amplitude. $`A`$ obeys the following equation

$$A_t=A+A_{xx}|A|^2A,$$
(64)
where coefficients can always be set to unity by an appropriate rescaling. This equation is universal in the sense that its form depends only on translational and rotational symmetries. It is usually referred to as the amplitude equation or the Ginzburg-Landau equation cross93 . The amplitude equation describes the nature of the bifurcation from the structureless state as well as instabilities of the ordered pattern with respect to wavelength modulations (the Eckhaus instability cross93 ).

A description in terms of a slowly varying amplitude makes a sense only if $`A`$ varies in a slow fashion in comparison with the pattern wavelength $`1/q_c`$. This requirement fails if $`q_c0`$. In that case a separation between a fast and slow variation is illegitimate. One has thus to resort to a (singular) expansion à la Sivashinsky sivashinsky77 . Because $`q_c`$ is small, it provides an appropriate parameter of expansion. This is also known as the long-wavelength limit. In that limit the field $`h`$ obeys generically the Kuramoto-Sivashinsky equation (or the damped form with $`\lambda `$ the damping rate), which is written in the canonical form,

$$h_t=\lambda hh_{xx}h_{xxxx}+h_x^2.$$
(65)
Again this equation is universal misbah94a ; remark  in the limit $`q_c0`$. Because of the smallness of $`q_c`$ (long-wavelength) any nonlocal effect present in the original constitutive equations disappears, and one recovers again a local equation. This is true only because a length scale is present and serves as a cut-off. It is akin to the Debye length which makes the nonlocal response function of the Coulomb gas local. There is, however, a myriad of situations where long-range interactions play a decisive role at all scales, and where it is hardly believable that the above equation (or a similar equation) should arise even in the long-wavelength limit. Typical situations with long range interactions are electrostatic and magnetic systems, problems with elastic fields, and so on. These systems are devoid of an intrinsic length-scale and it is highly desirable to derive the appropriate generic form of the corresponding evolution equation in the long-wavelength limit. A derivation on the basis of a gradient expansion of an appropriately chosen general nonlinear operator in the spirit of misbah94a  will be given elsewhere kassner02 .

Here we argue that we can essentially read off the form of that equation from our result (62). To this end, we simply expand to second order in $`\zeta _x`$ and drop all higher-order terms. As noted before, the kept terms are exact to that order. The “generic” evolution equation with long range interaction obtained this way reads

$$\zeta _t=\alpha \zeta +\zeta _{xx}\frac{2}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _x(x^{})}{x^{}x}dx^{}2\zeta _x^2\frac{1}{2}\left(\frac{2}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _x(x^{})}{x^{}x}dx^{}\right)^2.$$
(66)
Comparing this equation with the KS one (65), we see that it is not exactly a generalization. Since the nonlocal term drives the instability, the diffusive term need not be negative and hence there is no need for a fourth-order term. Rather, Eq. (66) is complementary to the Kuramoto-Sivashinsky result. Of course, the analog of Eq. (66) with a negative sign of $`\zeta _{xx}`$ may also appear and then the dispersion relation would need a different ultraviolet stabilization. If this term were fourth order in $`q`$, we would obtain a linear term $`\zeta _{xxxx}`$ and an equation that is a common generalization of the KS equation and Eq. (66).

Closer inspection shows that this equation can be only a particular case of the generic equation. For in general the amplitude equation will have coefficients in front of each of its six terms. Four of these can be made equal to one via division by a common prefactor and rescaling of space, time and the amplitude $`\zeta `$ itself. So our final equation should contain two independent nondimensional parameters, but it has only one, $`\alpha `$. We suspect that this is due to a hidden symmetry, and a natural candidate for this symmetry is the variational nature of the basic evolution dynamics (1,3). In the fully generic equation, the additional parameter would be the ratio of the coefficients of the two nonlinear terms. It is therefore easy to conjecture the form of this equation:

$$\zeta _t=\alpha \zeta +\zeta _{xx}\frac{2}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _x(x^{})}{x^{}x}dx^{}2\zeta _x^2\frac{\beta }{2}\left(\frac{2}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _x(x^{})}{x^{}x}dx^{}\right)^2,$$
(67)
where we have introduced a second parameter $`\beta `$. The corresponding universal equation for conservative dynamics would be obtained by preceding the complete right-hand side of this equation by the operator $`_x^2`$.

We have simulated Eq. (66) and find that it still describes coarsening, this time without saturation. An example of a simulation is displayed in Fig. 4. We see a single groove or “finger” survive after a long time.

On physical grounds, an non-saturating amplitude is to be expected in some situations involving the Grinfeld instability kassner01 . Nevertheless, even in such a case one should a priori not believe that a long-wave equation predicting an ever-increasing amplitude correctly describes the long-time asymptotics, because the interface eventually leaves the domain of validity of the long-wave assumption. This happens when the amplitude becomes larger than the wavelength. Therefore, an equation such as (67) should be expected to describe some intermediate-time asymptotics at best.

Things are different for other universal equations such as the Ginzburg-Landau equation (64) and the (damped) Kuramoto-Sivashinsky equation (65), which usually do work in the long-time limit. But this is due to the fact that in these cases the final (average) amplitude does saturate and is predicted to be small. However, if the cubic term in (64) is positive, which necessitates a fifth-order describtion to obtain saturation of the amplitude, the final asymptotics cannot ordinarily be trusted either, because then the amplitude may not be small contrary to the original assumptions.

On the other hand, while any amplitude description leading to non-saturating amplitudes may seem self-defeating in this sense, there appear to be systems, where the domain of validity of the amplitude equation extends beyond that of the assumptions made in its derivation. Recently, a singular long-wave equation was derived for step-flow growth without desorption pierrelouis98 . It predicts an amplitude growing as the square root of time, i.e., without bounds. Currently available evidence on the basis of numerical solutions of the exact model equations seems to corroborate that prediction. Hence, it is imperative to compare the long-time limit of the amplitude description to whatever information is available about the long-time behavior of the full system, in order to assess its validity.

Equation (67) should be thought of as being the first in a hierarchy of generic equations describing the initial and intermediate stages of evolution of a number of different systems. It is relevant to systems the linear dispersion relation of which is quadratic. For other systems, such as the Saffman-Taylor system, where the dispersion relation is cubic, we will have a different linear term in the equation (the Laplacian of a Hilbert transform), and it is an interesting question whether the appearing nonlinearities will remain the same. Following the same strategy for the Saffman-Taylor problem, we arrive at the following equation

$$\zeta _t=\frac{V}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _x(x^{})}{x^{}x}dx^{}\frac{\gamma B}{\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\frac{\zeta _{xxx}(x^{})}{x^{}x}dx^{}+V\zeta _x^2+\gamma B\zeta _x\zeta _{xxx},$$
(68)
where $`V`$ is the flat-interface velocity, $`\gamma `$ the surface tension, and $`B=b^2/12\mu `$ with $`b`$ the thickness of the Hele Shaw cell and $`\mu `$ the fluid viscosity. Of course the linear terms are always quite trivial, but what is noteworthy is that $`\beta =0`$, i.e., the nonlinear terms are local, like in the KS equation. (Under the usual assumptions for long-wave equations, the second local nonlinearity is small in comparison with the KS one.) We call this limit the Sivashinsky limit where nonlocality is present in the linear terms only, as was derived originally by Sivashinsky for flame propagation sivashinsky77 . This contrasts with our equation (67) which involves in the long-wavelength limit a nonlinear and nonlocal term. Therefore, this equation should introduce a new universality class with long range interactions.

## VI conclusions

To summarize, our method of derivation is applicable whenever the general solution of the bulk equations can be found by a transformation method. Usually this means they must be linear (with constant coefficients). A second condition, requiring that the dispersion relation can be transformed to long-wave form, should be always satisfiable. The generalization to three dimensions is straightforward, in principle. Our approach opens up the road to a realm of amplitude equations constituting a hitherto unknown type, the investigation of which seems worthwhile both from the mathematical and physical points of view; a new line of research should be stimulated.

The fact that the inner expansion produces only polynomials in $`Z`$ (and hence the outer solution is identical to the uniform approximation Bender78 ) strongly suggests that the derivation may even be achieved without the detour via a long-wave approach, using a regular perturbation expansion of the outer equations. This would seem to indicate that we are not really dealing with an equation that is valid only in the long-wave limit but rather with an ordinary (but strongly nonlinear) amplitude equation valid at all wavelengths for small enough amplitudes, thus a posteriori justifying the formal device of rescaling to large wavelengths. On the other hand, the first of our approximations corresponds to setting $`\left|q\right|\zeta `$ equal to zero in the exponents of (27) and (28), which, after taking into account all rescalings, means that the amplitude is small compared with the wavelength.

Finally, the amplitude equation is much faster (and easier) to handle numerically than previous schemes. Like the KS equation, for example, its leading behavior described by Eq. (67) should be generic for systems undergoing long-wavelength instabilities with long-range interactions.

Acknowledgment We are grateful to Peter Kohlert for help with the computer algebra.


format: text

text:
# Focusing Vacuum Fluctuations II

## I Introduction

In a previous paper FS00 , henceforth I, we developed a geometric optics approach to the quantization of scalar and electromagnetic fields near the focus of a parabolic mirror. We found that there can be enhanced fluctuations near the focus in the sense that mean squared field quantities scale as an inverse power of the distance from the focus, rather than an inverse power of the distance from the mirror. These enhanced fluctuations were found to arise from an interference term between different reflected rays. In the present paper we extend the previous treatment. In I, only points on the symmetry axis were considered. Here we are able to treat points in an arbitrary direction from the focal line of a parabolic cylinder. We give more detailed numerical results which provide a fuller picture of the phenomenon of focusing of vacuum fluctuations. We also correct some erroneous results in Sect. V of I.

The outline of the present paper is as follows: In Sect. II, we review and extend some of the formalism used to calculate the mean squared scalar and electric fields, $`\phi ^2`$ and $`𝐄^2`$, near the focus. In Sect. III we derive some geometric expressions which are needed to study fluctuations at points off of the symmetry axis. In Sect. IV, the evaluation of integrals with singular integrands is revisited. Two different, but equivalent, approaches are discussed. A particular case where the integrals can be performed analytically is examined in Sect. V. More generally, it is necessary to calculate $`\phi ^2`$ and $`𝐄^2`$ numerically. A procedure for doing so is outlined in Sect. VI. Some detailed numerical results are also presented there. The limits of validity of our model and results will be examined in Sect. VII. This discussion will draw on some results on diffraction obtained in the Appendix. The experimental testability of our conclusions will be discussed in Sect. VIII. Finally, the results of the paper will be summarized in Sect. IX.

Units in which $`\mathrm{}=c=1`$ will be used throughout this paper. Electromagnetic quantities will be in Lorentz-Heaviside units.

## II Basic Formalism

Here we will briefly review the geometric optics approach developed in I. The basic assumption is that we may use a ray tracing method to determine the functional form of the high frequency modes, which will in turn give the dominant contribution to the expectation values of squared field operators. We start with a basis of plane wave modes. The incident wave, for a scalar field, may be taken to be

$$f_𝐤=\frac{1}{\sqrt{2\omega V}}\mathrm{e}^{i(𝐤𝐱\omega t)},$$
(1)
with box normalization in a volume $`V`$. In the presence of a boundary, this is replaced by the sum of incident and reflected waves,

$$F_𝐤=f_𝐤+\underset{i}{}f_𝐤^{(i)},$$
(2)
where the $`f_𝐤^{(i)}`$ are the various reflected waves. One could also in principle adopt a wavepacket basis, in which $`F_𝐤`$ is replaced by a localized wavepacket. Because the time evolution preserves the Klein-Gordon norm, if the various modes are orthonormal in the past, they will remain so after reflection from the mirror. Thus we can view Eq. (2) as the limit of a set of orthonormal wavepacket modes in which the modes become sharply peaked in frequency and hence delocalized.

It was shown in I that the renormalized expectation value of the squared scalar field is given by a sum of interference terms:

$$\phi ^2=\underset{𝐤}{}\left[\underset{i}{}(f_𝐤^{}f_𝐤^{(i)}+f_𝐤f_{}^{(i)}{}_{𝐤}{}^{})+\underset{ij}{}f_𝐤^{(i)}f_{}^{(j)}{}_{𝐤}{}^{}\right].$$
(3)
This renormalized expectation value is defined as a difference in the mean value of $`\phi ^2`$ with and without the mirror, and hence will vanish at large distances from the mirror. The various interference terms in the above expression yield contributions to $`\phi ^2`$ which are proportional to the inverse square of the appropriate path difference. Thus in the vicinity of the focus, the interference terms between different reflected rays will dominate over that between the incident and a reflected ray. In the present paper, we will consider cases with no more than two reflected rays, and write

$$\phi ^22\mathrm{Re}\underset{𝐤}{}f_𝐤^{(1)}f_{}^{(2)}{}_{𝐤}{}^{}.$$
(4)
In the case of a parabolic cylinder, this may be expressed as

$$\phi ^2=\frac{1}{3\pi ^2}\frac{d\theta ^{}}{(\mathrm{\Delta }\mathrm{})^2}.$$
(5)
Here $`\mathrm{\Delta }\mathrm{}`$ is the path difference for the two reflected rays, and the integration is over the reflection angle of one of the rays. The range is chosen so that each pair of reflected rays is counted once. The corresponding expression for the mean squared electric field near the focus of a parabolic cylinder is found in I to be

$$𝐄^2=\frac{8}{5\pi ^3}\frac{d\theta ^{}}{(\mathrm{\Delta }\mathrm{})^4}.$$
(6)
## III Optics of Parabolic Mirrors

In this section, we wish to generalize some of the results of I concerning the incident and reflected rays in the presence of a parabolic mirror. Consider the geometry illustrated in Fig. 1. An incident ray at an angle of $`\theta `$ to the symmetry axis is reflected from the point $`(x_i,y_i)`$, and then reaches the point $`P`$ at an angle of $`\theta ^{}`$. We first need to find the relation between $`\theta `$ and $`\theta ^{}`$. Note that the reflected ray crosses the symmetry axis at a distance $`c`$ from the focus. It was shown in I that

$$\theta =\frac{c\mathrm{sin}^3\theta ^{}}{b(1\mathrm{cos}\theta ^{})}.$$
(7)
However, from the law of sines, we have that

$$c\mathrm{sin}\theta ^{}=a\mathrm{sin}(\theta ^{}\gamma ).$$
(8)
Hence we now have that

$$\theta =\frac{a}{b}f(\theta ^{}),$$
(9)
where

$$f(\theta ^{})=\frac{\mathrm{sin}^3\theta ^{}\mathrm{sin}(\theta ^{}\gamma )}{(1\mathrm{cos}\theta ^{})}.$$
(10)
There will be multiply reflected rays whenever different values of $`\theta ^{}`$ are associated with the same value of $`f`$. The function $`f(\theta ^{})`$ is plotted in Fig. 2 for various values of $`\gamma `$. We can see from these plots that in general there can be up to four reflected angles $`\theta ^{}`$ for a given incident angle $`\theta `$. However, if the mirror size $`\theta _0`$ is restricted to be less than $`2\pi /3`$, then there will never be more than two values of $`\theta ^{}`$ for a given $`\theta `$. Throughout this paper, we will assume $`\theta _0<2\pi /3`$, and hence have at most two reflected rays for a given incident ray. The two reflected rays will occur at $`\theta ^{}=\alpha `$ and $`\theta ^{}=\beta `$, where

$$f(\alpha )=f(\beta ).$$
(11)
Our next task is to calculate the difference in path length, $`\mathrm{\Delta }\mathrm{}`$, for these two reflected rays. Again, this is a generalization of a calculation given in I. Consider the situation illustrated in Fig. 3, where a reflected ray with angle $`\theta ^{}`$ reaches the point $`P`$ after reflecting from the point $`(x_i,y_i)`$ on the mirror. Let $`s_1`$ be the distance traveled after reflection, and $`s_2`$ be the distance traveled between when the ray crosses the line $`x=x_0`$ and when it reaches the mirror. Note that

$$s_2=(x_ix_0)\mathrm{sec}\theta x_ix_0+O\left(\frac{a^2}{b^2}\right).$$
(12)
Even if $`x_i<0`$, we can choose $`x_0`$ to be such that $`s_2>0`$. Next note that

$$s_1=(y_ia\mathrm{sin}\gamma )\mathrm{csc}\theta ^{}.$$
(13)
The reflection point $`(x_i,y_i)`$ is the intersection of the line

$$y=\mathrm{tan}\theta ^{}(xa\mathrm{cos}\gamma )+a\mathrm{sin}\gamma $$
(14)
with the parabola, given by

$$x=\frac{b^2y^2}{b}.$$
(15)
These relations lead to

$$y_i=b\mathrm{cot}\theta ^{}\left[\mathrm{sec}\theta ^{}\sqrt{12\left(\frac{a}{b}\right)\mathrm{sin}^2\theta ^{}(\mathrm{cos}\gamma \mathrm{sin}\gamma \mathrm{cot}\theta ^{})}1\right].$$
(16)
We may now expand this expression to first order in $`a/b`$ and combine it with our previous expressions to show that, to first order,

$$\mathrm{}=s_1+s_2=bx_0a(\mathrm{cos}\gamma \mathrm{cos}\theta ^{}\mathrm{sin}\gamma \mathrm{sin}\theta ^{}).$$
(17)
Thus the magnitude of the path length difference for $`\alpha `$ and $`\beta `$, two different values of $`\theta ^{}`$, is

$$\mathrm{\Delta }\mathrm{}=a|\mathrm{cos}\gamma (\mathrm{cos}\alpha \mathrm{cos}\beta )+\mathrm{sin}\gamma (\mathrm{sin}\alpha \mathrm{sin}\beta )|.$$
(18)
In the limit that $`\gamma =0`$, we obtain the result for $`\mathrm{\Delta }\mathrm{}`$ used in I.

## IV Evaluation of Singular Integrals

We can express Eq. (5) as

$$\phi ^2=\frac{1}{6\pi ^2a^2}\frac{d\alpha }{h^2},$$
(19)
and Eq. (6) as

$$𝐄^2=\frac{4}{5\pi ^3a^4}\frac{d\alpha }{h^4}.$$
(20)
Here

$$h=h(\alpha ,\beta )=|\mathrm{cos}\gamma (\mathrm{cos}\alpha \mathrm{cos}\beta )+\mathrm{sin}\gamma (\mathrm{sin}\alpha \mathrm{sin}\beta )|,$$
(21)
and the integrations run over the full range for which there are two reflected rays, thus counting each pair twice.

The integrands of these integrals contain singularities within the range of integration. These singularities occur at a critical angle, $`\theta _c`$, at which $`h`$ vanishes linearly. Thus the $`\phi ^2`$ integral contains a $`(\alpha \theta _c)^2`$ singularity, and the $`𝐄^2`$ contain a $`(\alpha \theta _c)^4`$ singularity. These singularities are presumably artifacts of assuming perfect mirrors with sharp boundaries as will be discussed in Sect. VII. The singularity occurs when both $`\alpha `$ and $`\beta `$ are approaching $`\theta _c`$ from opposite directions. The critical angles $`\theta _c`$ are just the extrema of the function $`f(\theta ^{})`$.

Integrals with such singular integrands can be defined by a generalization of the principal value prescription. This generalization involves an integration by parts to recast the original integral as one containing a less singular integrand, plus surface terms. For example, we may use

$$\frac{1}{x^2}=\frac{1}{2}\frac{d^2}{dx^2}\mathrm{ln}x^2$$
(22)
to write

$$_a^b𝑑x\frac{f(x)}{x^2}=\frac{1}{2}_a^b𝑑xf^{\prime \prime }(x)\mathrm{ln}x^2+\left[\frac{f(x)}{x}+\frac{1}{2}f^{}(x)\mathrm{ln}x^2\right]_a^b.$$
(23)
Similarly,

$$\frac{1}{x^4}=\frac{1}{12}\frac{d^4}{dx^4}\mathrm{ln}x^2,$$
(24)
leads to

$$_a^b𝑑x\frac{f(x)}{x^4}=\left[\frac{f(x)}{3x^3}+\frac{f^{}(x)}{6x^2}+\frac{f^{\prime \prime }(x)}{6x}\frac{1}{12}f^{\prime \prime \prime }(x)\mathrm{ln}x^2\right]_a^b\frac{1}{12}_a^b𝑑x\mathrm{ln}x^2\frac{d^4f(x)}{dx^4}.$$
(25)
If $`b>0`$ and $`a<0`$, there were singularities in the original integrals which are replaced by integrable, logarithmic singularities. In all cases, the surface terms are evaluated away from the singularity and are hence finite. Note that we could have added additional polynomial terms on the righthand sides of both Eqs. (22) and (24). For example, we could replace $`\mathrm{ln}x^2`$ by $`\mathrm{ln}x^2+c_1x+c_0`$ in Eq. (22). However, the arbitrary constants $`c_0`$ and $`c_1`$ will cancel out between the two terms on the righthand side of Eq. (23) and hence can be ignored. A similar cancellation occurs if we add additional terms into Eq. (24).

If $`f`$ and a sufficient number of its derivatives vanish at the endpoints, then the surface terms vanish. In I, it was incorrectly argued that the surface terms which arise in the present problem can be ignored. It is true that if the reflectivity of the mirror falls smoothly to zero, the surface terms vanish. However, in this case it would be necessary to integrate explicitly the function which is falling smoothly to zero. If the reflectivity falls very rapidly at the edge of the mirror, its derivatives will be large, and effectively reproduce the surface terms. Thus, the explicit results given in Sect. 5 of I are incorrect. In Sects. V and VI of the present paper, we will give new results which replace and generalize the older results.

There is an alternative way to implement the above integration by parts prescription. This is simply to evaluate the integral as an indefinite integral, and evaluate the result at the endpoints, ignoring the singularity within the integration range. Even if it is not possible to find a closed form expression for the indefinite integral, one can expand the integrand in a Laurent series about $`x=0`$ and integrate term by term, using the relation

$$_a^b\frac{1}{x^n}𝑑x=\frac{1}{n+1}\left(\frac{1}{a^{n+1}}\frac{1}{b^{n+1}}\right).$$
(26)
One can verify this relation for the case $`a<0<b`$ using the above integration by parts method.

The integrands in Eqs. (19) and (20) arise from the integrals

$$\frac{1}{x^2}=_0^{\mathrm{}}𝑑\omega \omega \mathrm{cos}\omega x$$
(27)
and

$$\frac{1}{x^4}=\frac{1}{6}_0^{\mathrm{}}𝑑\omega ^3\omega \mathrm{cos}\omega x,$$
(28)
respectively. The integrals are understood to be performed with the aid of a convergence factor, such as $`\mathrm{e}^{\lambda \omega }`$, with the limit $`\lambda 0`$ taken after integration. The singularities at $`x=0`$ are due to the contributions of arbitrarily large values of $`\omega `$. If there were to be a physical cutoff at high frequencies, then these singulariries would disappear. However, the integrals containing the singular integrands are independent of the cutoff in the limit that it occurs at sufficiently high frequencies. As an example, consider the function

$$g_2(x,\lambda )=_0^{\mathrm{}}𝑑\omega \omega \mathrm{e}^{\lambda \omega }\mathrm{cos}\omega x=\frac{\lambda ^2x^2}{(\lambda ^2+x^2)^2}.$$
(29)
As $`\lambda 0`$, $`g_2(x,\lambda )1/x^2`$, but for $`\lambda 0`$, $`g_2(x,\lambda )`$ is finite for all $`x`$. Furthermore, if we integrate $`g_2(x,\lambda )`$ on $`x`$ through $`x=0`$, and then take the $`\lambda 0`$ limit, the result will be finite and the same as that obtained by the above formal procedures:

$$_{x_0}^{x_0}g_2(x,\lambda )𝑑x=\frac{2x_0}{\lambda ^2+x_0^2}\frac{2}{x_0}.$$
(30)
Similarly, if we define

$$g_4(x,\lambda )=_0^{\mathrm{}}𝑑\omega \omega ^3\mathrm{e}^{\lambda \omega }\mathrm{cos}\omega x=\frac{6(x^22\lambda x\lambda ^2)(x^2+2\lambda x\lambda ^2)}{(\lambda ^2+x^2)^4},$$
(31)
then $`g_4(x,\lambda )6/x^4`$ as $`\lambda 0`$. However,

$$_{x_0}^{x_0}g_4(x,\lambda )𝑑x=\frac{4x_0(x_0^23\lambda ^2)}{\lambda ^2+x_0^3}\frac{4}{x_0^3}.$$
(32)
Thus we can understand the physical meaning of the formal integration procedures discussed above as follows: they provide a shortcut method to obtain the results one would obtain by inserting a frequency cutoff in the $`\omega `$ integrals, and then removing this cutoff after integration over $`x`$.

One peculiar feature of our results will be that both $`\phi ^2`$ and $`𝐄^2`$ diverge for particular values of $`\gamma `$. This arises when $`\theta _c`$, the extremum of $`f(\theta ^{})`$, approaches the edge of the mirror, $`\theta ^{}=\pm \theta _0`$. The mathematical reason for the divergence is that one limit of integration is approaching a point at which $`h=0`$ and the integrand is singular. This corresponds to letting either $`a0`$ with $`b`$ fixed, or $`b0`$ with $`a`$ fixed, in Eq. (26). The physical origin of this singular behavior is that we have made two unrealistic assumptions. The mirror is assumed to be perfectly reflecting at all frequencies and to have sharp edges at $`\theta ^{}=\pm \theta _0`$. This will be discussed further in Sect. VII.

## V Exact Results for $`\gamma =\frac{\pi }{2}`$

In general, the integrals for $`\phi ^2`$ and $`𝐄^2`$, Eqs. (5) and (6), respectively, can only be evaluated numerically. This is in part because the second reflection angle $`\beta `$ is implicitly given as a function of the first angle, $`\alpha `$. There is one case in which this relation may be written down in closed form. This is when $`\gamma =\frac{\pi }{2}`$, so the function $`f(\theta ^{})`$ is symmetrical about the origin. In this case, we have

$$\beta =\alpha ,$$
(33)
and we can write

$$\phi ^2=\frac{1}{24\pi ^3a^2}_{\theta _0}^{\theta _0}\frac{d\alpha }{\mathrm{sin}^2\alpha }=\frac{\mathrm{cot}\theta _0}{12\pi ^3a^2}.$$
(34)
Similarly,

$$𝐄^2=\frac{1}{20\pi ^3a^4}_{\theta _0}^{\theta _0}\frac{d\alpha }{\mathrm{sin}^4\alpha }=\frac{\mathrm{cos}\theta _0(32\mathrm{cos}^2\theta _0)}{30\pi ^3a^4\mathrm{sin}^3\theta _0}.$$
(35)
These expressions are plotted in Fig. 4. Note that both $`\phi ^2`$ and $`𝐄^2`$ can have either sign:

$`\phi ^2>0\mathrm{and}𝐄^2<0\mathrm{for}0<\theta _0<{\displaystyle \frac{\pi }{2}},`$
$`\phi ^2<0\mathrm{and}𝐄^2>0\mathrm{for}{\displaystyle \frac{\pi }{2}}<\theta _0<{\displaystyle \frac{2\pi }{3}}.`$ (36)
Both quantities seem to diverge in the limit that $`\theta _00`$

$`\phi ^2{\displaystyle \frac{1}{12\pi ^3a^2\theta _0}},`$
$`𝐄^2{\displaystyle \frac{1}{30\pi ^3a^4\theta _0^3}}.`$ (37)
This divergence is a special case of the singularities arising when the edge of the mirror approaches an extremum of $`f(\theta ^{})`$.

We can go further and find the derivatives of both $`\phi ^2`$ and $`𝐄^2`$ with respect to $`\gamma `$, at $`\gamma =\frac{\pi }{2}`$. First, we need to find $`\beta (\alpha )`$ to first order in $`\gamma \frac{\pi }{2}`$. We may expand the relation $`f(\beta )=f(\alpha )`$ to this order and show that

$$\beta =\alpha \frac{2\mathrm{sin}^2\alpha (\gamma \frac{\pi }{2})}{(\mathrm{cos}\alpha 1)(2\mathrm{cos}\alpha +1)}+\mathrm{}.$$
(38)
We can next expand $`1/h^2`$ and $`1/h^4`$ to first order in $`\gamma \frac{\pi }{2}`$. The first order corrections are odd functions of $`\alpha `$ whose explicit forms we will not need.

The crucial effect of varying $`\gamma `$ slightly away from $`\pi /2`$ is to change the range of integration. First consider the case where $`\gamma <\frac{\pi }{2}`$. The integration range now becomes

$$\theta _0<\alpha <\theta _0\delta ,$$
(39)
where

$$\delta =\frac{2\mathrm{sin}^2\theta _0(\frac{\pi }{2}\gamma )}{(1\mathrm{cos}\theta _0)(2\mathrm{cos}\theta _0+1)}.$$
(40)
To first order in $`\gamma \frac{\pi }{2}`$, we can write

$`\phi ^2`$ $`=`$ $`{\displaystyle \frac{1}{6\pi ^3a^2}}{\displaystyle _{\theta _0}^{\theta _0\delta }}{\displaystyle \frac{d\alpha }{h^2}}`$ (41)
$`=`$ $`{\displaystyle \frac{1}{6\pi ^3a^2}}\left[{\displaystyle _{\theta _0}^{\theta _0}}𝑑\alpha \left({\displaystyle \frac{1}{h^2}}\right)_0+{\displaystyle _{\theta _0}^{\theta _0}}𝑑\alpha \left({\displaystyle \frac{1}{h^2}}\right)_1{\displaystyle _{\theta _0\delta }^{\theta _0}}𝑑\alpha \left({\displaystyle \frac{1}{h^2}}\right)_0\mathrm{}\right],`$
where the subscripts on $`1/h^2`$ denote the order in an expansion in powers of $`\gamma \frac{\pi }{2}`$. The first term on the right hand side of Eq. (41) is just the zeroth order part given by Eq. (34). The next term vanishes because $`(1/h^2)_1`$ is an odd function. The final term may be approximated using

$$_{\theta _0\delta }^{\theta _0}𝑑\alpha \left(\frac{1}{h^2}\right)_0\delta \left[\left(\frac{1}{h^2}\right)_0\right]_{\alpha =\theta _0}.$$
(42)
This leads to

$$\phi ^2=\frac{1}{12\pi ^3a^2}\left[\frac{1}{\mathrm{tan}\theta _0}\frac{\frac{\pi }{2}\gamma }{(1\mathrm{cos}\theta _0)(2\mathrm{cos}\theta _0+1)}+\mathrm{}\right],\gamma <\frac{\pi }{2}.$$
(43)
We may now repeat this procedure for $`\gamma >\frac{\pi }{2}`$. In this case, the range of integration is

$$\theta _0\delta <\alpha <\theta _0,$$
(44)
and $`\delta <0`$, so the range of integration has again decreased. This causes the first order change in $`\phi ^2`$ to have the same magnitude but the opposite sign from the previous case. Thus in all cases, we can write

$$\phi ^2=\frac{1}{12\pi ^3a^2}\left[\frac{1}{\mathrm{tan}\theta _0}\frac{|\gamma \frac{\pi }{2}|}{(1\mathrm{cos}\theta _0)(2\mathrm{cos}\theta _0+1)}+\mathrm{}\right].$$
(45)
A similar analysis may be applied to $`𝐄^2`$ with the result

$$𝐄^2=\frac{4}{5\pi ^3a^4}\left[\frac{\mathrm{cos}\theta _0(32\mathrm{cos}^2\theta _0)}{24\mathrm{sin}^3\theta _0}\frac{|\gamma \frac{\pi }{2}|}{8\mathrm{sin}^2\theta _0(1\mathrm{cos}\theta _0)(2\mathrm{cos}\theta _0+1)}+\mathrm{}\right].$$
(46)
Thus both $`\phi ^2`$ and $`𝐄^2`$ have cusps at $`\gamma =\frac{\pi }{2}`$. The nonanalytic behavior is due to the fact that the range of integration decreases whenever $`\gamma `$ moves away from $`\frac{\pi }{2}`$ in either direction. As will be discussed in Sect. VII, the cusp is presumably an artifact of a sharp edge approximation, and should be smoothed out in a more exact treatment.

## VI Numerical Procedures and Results

Apart from the special cases discussed in the previous section, it is necessary to evaluate $`\phi ^2`$ and $`𝐄^2`$ numerically. The first step is to find the second reflection angle $`\beta `$ as a function of the first reflection angle $`\alpha `$. This involves a straightforward numerical solution of the equation

$$f(\beta )=f(\alpha ).$$
(47)
Because we assume a restriction on the angle size of the mirror,

$$\theta _0<\frac{2\pi }{3},$$
(48)
there will always be either one or no roots for $`\beta `$. Within the geometric optics approximation that we use, integrands are assumed to vanish in regions where there are no roots.

As was discussed in Sect. IV, there are at least two methods that may be used for explicit evaluation of the integrals which appear in the geometric optics expressions for $`\phi ^2`$ and $`𝐄^2`$. The first is an integration by parts, which replaces the singularity in Eqs. (19) and (20) by a logarithmic singularity. Here we will outline how this may be done explicitly. Consider first Eq. (19), which may be expressed as

$$\phi ^2=\frac{1}{6\pi ^2a^2}\frac{d\alpha }{[h(\alpha ,\beta (\alpha ))]^2}.$$
(49)
Suppose that we are interested in integrating over the range $`\alpha _1<\alpha <\alpha _2`$ and that this range contains one point $`\alpha =\alpha _0`$ at which $`h=0`$. We can write

$$\phi ^2=\frac{1}{12\pi ^2a^2}_{\alpha _1}^{\alpha _2}d\alpha \frac{(\alpha \alpha _0)^2}{[h(\alpha ,y(\alpha ))]^2}\frac{d^2}{d\alpha ^2}\mathrm{log}(\alpha \alpha _0)^2.$$
(50)
The quantity $`(\alpha \alpha _0)^2/[h(\alpha ,y(\alpha ))]^2`$ is finite at $`\alpha =\alpha _0`$. This expression may now be integrated by parts, as in Eq. (23). The derivatives of $`h`$ which arise are computed using Eq. (21) and the relation

$$\frac{d\beta }{d\alpha }=\frac{f^{}(\alpha )}{f^{}(\beta )}.$$
(51)
An analogous procedure can be applied to $`𝐄^2`$. The detailed expressions for $`\phi ^2`$ and $`𝐄^2`$ which result from this procedure are rather complicated, and will not written down explicitly.

The second method which may be employed is a variant of the direct integration illustrated in Eq. (26). The actual integrands in Eqs. (19) and (20) are too complex to integrate in closed form. However, we can use numerical integration of $`1/h^2`$ and $`1/h^4`$ in regions away from zeros of $`h`$ and direct integration of a series expansion in the neighborhood of a zero. The first step in the generation of the series expansion is to expand $`\beta (\alpha )`$ around $`\alpha =\alpha _0`$,

$$\beta =2\alpha _0\alpha +a_2(\alpha \alpha _0)^2+a_3(\alpha \alpha _0)^3+\mathrm{}.$$
(52)
The coefficients are found by inserting this expansion into Eq. (47) and then expanding both sides of the resulting expression in powers of $`\alpha \alpha _0`$. A few of the leading coefficients are

$`a_2`$ $`=`$ $`{\displaystyle \frac{f^{\prime \prime \prime }}{3f^{\prime \prime }}},`$
$`a_3`$ $`=`$ $`{\displaystyle \frac{(f^{\prime \prime \prime })^2}{9(f^{\prime \prime })^2}},`$
$`a_4`$ $`=`$ $`{\displaystyle \frac{9(f^{\prime \prime })^2(f^{(5)})^230f^{\prime \prime }f^{\prime \prime \prime }f^{(4)}+40(f^{\prime \prime \prime })^2}{540(f^{\prime \prime })^3}},`$ (53)
where the derivatives of $`f`$ are evaluated at $`\alpha =\alpha _0`$. The expansion for $`y(\alpha )`$ is then used to generate analogous expansions for $`1/h^2`$ and $`1/h^4`$, which are in turn integrated over the interval $`\alpha _0\xi _0<\alpha <\alpha _0+\xi _0`$ using Eq. (26). The result is combined with the direct numerical integration outside of this interval. Here $`\xi _0`$ is an arbitrary small positive number. One of the tests of the numerical procedure is the independence of the results upon the choice of $`\xi _0`$.

We have developed numerical routines based upon both of the above procedures. The plots which are given below were created using a routine based upon the second method, with an expansion of $`y`$ to sixth order in $`\alpha \alpha _0`$. The particular values of $`\xi _0`$ were in the range $`0.1\xi _00.3`$. In this range, the routine is relatively stable and insensitive to $`\xi _0`$. Smaller values of $`\xi _0`$ can lead to instability, as the contributions from $`|xx_0|>\xi _0`$ and from $`|xx_0|<\xi _0`$ are both large in magnitude and tending to cancel each other.

Numerical results for $`\phi ^2`$ and $`𝐄^2`$ for various values of $`\theta _0`$ are shown in Figs. 5-12 as functions of $`\gamma `$. In all cases, there are values of $`\gamma `$ at which $`\phi ^2`$ or $`𝐄^2`$ are singular in the model of a perfectly reflecting mirror with sharp edges. These points occur when an extremum of the function $`f(\theta )`$ sits at an edge of the mirror, $`\theta ^{}=\pm \theta _0`$. The values and slopes of $`\phi ^2`$ and $`𝐄^2`$ at $`\gamma =\frac{\pi }{2}`$ are in agreement with Eqs. (45) and (46). We can see from the graphs that for smaller mirrors, $`\theta _0<\frac{\pi }{2}`$, we have $`\phi ^2>0`$ and $`𝐄^2<0`$ everywhere. For larger mirrors, $`\theta _0>\frac{\pi }{2}`$, we have regions where $`\phi ^2>0`$ and other regions where it is negative, and similarly for $`𝐄^2`$. The signs of $`\phi ^2`$ and of $`𝐄^2`$ always seem to be opposite. The regions near $`\gamma =\frac{\pi }{2}`$ for larger mirrors where $`𝐄^2>0`$ are of special interest. These are regions where an atom will feel an attractive force toward the focus, and thus has the possibility of being trapped.

## VII Limits of Validity of the Results

In this section, we will discuss the likely ranges of validity of the model we have used to calculate $`\phi ^2`$ and $`𝐄^2`$. In particular, we have assumed a geometric optics approximation, and a mirror which is both perfectly reflecting and has sharp edges. Each of these assumptions will be examined critically.

### VII.1 Geometric Optics

The use of geometric optics amounts to ignoring diffraction effects, so we can gauge the accuracy of geometric optics by estimating the size of these effects. In the geometric optics approximation, the incident and reflected waves in Eq. (2) all have the same magnitude. Suppose that we now introduce a correction due to diffraction and write

$$F_𝐤=f_𝐤+\underset{i}{}f_𝐤^{(i)}+\mathrm{\Delta }F_𝐤.$$
(54)
In the case of a plane strip which has width $`2y_0`$ in one direction and is infinite in the other direction, the magnitude of the diffraction correction is estimated in the Appendix with the result $`|\mathrm{\Delta }F_𝐤/F_𝐤|\sqrt{\lambda /y_0}`$, where $`\lambda `$ is the wavelength of the incident wave. In our case, we have a parabolic cylinder characterized by the length scale $`b`$. We could imagine representing the parabolic cylinder by a set of strips which are infinite in the $`z`$-direction. The order of magnitude of the diffraction correction should be the same as for a single strip of width $`b`$. Thus we estimate that in the present case

$$\left|\frac{\mathrm{\Delta }F_𝐤}{F_𝐤}\right|\sqrt{\frac{\lambda }{b}}.$$
(55)
The wavelengths which give the dominant contribution to $`\phi ^2`$ and $`𝐄^2`$ near the focus are those of order $`a`$. Thus the interference term between $`\mathrm{\Delta }F_𝐤`$ and any of the other terms in Eq. (54) should yield contributions to $`\phi ^2`$ and $`𝐄^2`$ which are smaller than the dominant contribution by a factor of the order of $`\sqrt{a/b}`$. Thus we estimate the diffraction contribution to $`\phi ^2`$ to be of order

$$\mathrm{\Delta }\phi ^2\frac{1}{a^{\frac{3}{2}}b^{\frac{1}{2}}},$$
(56)
and that to $`𝐄^2`$ to be of order

$$\mathrm{\Delta }𝐄^2\frac{1}{a^{\frac{7}{2}}b^{\frac{1}{2}}}.$$
(57)
Note that the geometric optics results near the focus depend only upon $`a`$ and $`\theta _0`$, the angular size of the mirror, but not upon $`b`$, the linear dimension of the mirror. However, the diffraction correction decreases with increasing $`b`$ and hence can be made smaller for a larger, more distant mirror.

In classical optics, the effects of diffraction are normally of the order of the wavelength divided by the size of the object. This arises when one is looking at the power in the diffracted wave, which is given by the square of Eq. (55). In our case, there is a possible contribution from an interference term between the diffracted wave and the geometric optics contributions. If a more detailed calculation were to find that this term vanished, then our estimates for $`\mathrm{\Delta }\phi ^2`$ and $`\mathrm{\Delta }𝐄^2`$ would become smaller than Eqs. (56) and (57), respectively, by an additional factor of $`\sqrt{a/b}`$.

### VII.2 Edge Effects and Finite Reflectivity

Our model of the mirror is one in which it is not only perfectly reflecting, but also has a sharp edge at which the reflectivity falls discontinuously to zero. Both of these assumption are over simplifications to which the singular behavior of $`\phi ^2`$ and $`𝐄^2`$ may be attributed. Consider the sharp edge assumption. A more realistic model might have the reflectivity falling smoothly to zero over an angular interval of width $`\mathrm{\Delta }\theta `$ at the edges of the mirror. This would remove the singularities found above at specific values of $`\gamma `$. Recall that these singularities arise when an edge of the mirror, at $`\pm \theta _0`$, sits at an extremum of $`f(\theta ^{})`$. In the case of a smoothed edge, both $`\phi ^2`$ and $`𝐄^2`$ will be bounded for fixed $`a`$

$$|\phi ^2|<\frac{1}{a^2\mathrm{\Delta }\theta }$$
(58)
and

$$|𝐄^2|<\frac{1}{a^4(\mathrm{\Delta }\theta )^3}.$$
(59)
Smoothing of the edges of the mirror is also expected to remove the cusps at $`\gamma =\frac{\pi }{2}`$.

The smoothed edges do not, however, remove the singularities as $`a0`$. This singularity is presumably due to the assumption of perfect reflectivity at all wavelengths. A more realistic model would have the reflectivity go to zero at short wavelengths. Suppose that the mirror becomes transparent for wavelengths less than some minimum value, $`\lambda _m`$. Then we expect to find the bounds

$$|\phi ^2|<\frac{1}{\lambda _m^2}$$
(60)
and

$$|𝐄^2|<\frac{1}{\lambda _m^4}.$$
(61)
One reason for reduced reflectivity at short wavelengths is dispersion. The mirror can be regarded as close to perfectly reflecting only for wavelengths longer than about the plasma wavelength $`\lambda _P`$ of the metal in the mirror. Thus our results are only valid when $`a>\lambda _P`$. For aluminum, for example, $`\lambda _P84\mathrm{n}\mathrm{m}`$. However, it is unlikely that dispersion alone is capable of removing all short distance singularities. The reason for this is that dielectric functions approach unity as $`(\lambda /\lambda _P)^2`$ as $`\lambda 0`$. This is not fast enough to regulate integrals which diverge quartically as the short wavelength cutoff. In the case of a plane interface between vacuum and a dispersive material, quantities such as $`𝐄^2`$ still diverge as the interface is approached, although less rapidly than in the case of a perfect mirror SF . Thus it seem that other effects, such as surface roughness, the breakdown of the continuum description at the atomic level, or quantum uncertainty in the location of the interface FS98  are needed to produce finite values of $`𝐄^2`$. If any effect does produce a sufficiently sharp cutoff at short wavelengths below about $`\lambda _m`$, then we would have the bounds in Eqs. (60) and (61) even without removing the assumption of a sharp edge to the mirror.

## VIII Observable Consequences

In I, several possible experimental tests of the enhanced vacuum fluctuation near the focus were discussed. Here these tests will be reviewed and further discussed. The most direct test would seem to be to measure the force on an atom or other polarizable particle. In a regime where the atom can be described by a static polarizability $`\alpha `$, the force is $`𝐅=V`$, where

$$V=\frac{1}{2}\alpha 𝐄^2.$$
(62)
In the vicinity of the focus, we have found that

$$𝐄^2=\frac{\mathrm{\Lambda }}{a^4},$$
(63)
where $`\mathrm{\Lambda }`$ is a dimensionless constant. For an atom in its ground state, Eq. (62) should be a good approximation when $`a`$ is larger than the wavelength associated with the transition to the first excited state.

One might try to measure the defection of atoms moving parallel to the focal line of a parabolic cylinder. The analogous experiment for a flat plate was performed by Sukenik et al Sukenik  and confirmed Casimir and Polder’s theoretical prediction CP . In the present case, the expected angular deflection is

$$\frac{\mathrm{\Delta }a}{a}=0.25\left(\frac{\mathrm{\Lambda }}{10^3}\right)\left(\frac{\alpha }{\alpha _{Na}}\right)\left(\frac{m_{Na}}{m}\right)\left(\frac{1\mu \mathrm{m}}{a}\right)^6\left(\frac{t}{10^3s}\right)^2.$$
(64)
Here $`m_{Na}=3.8\times 10^{23}\mathrm{gm}`$ and $`\alpha _{Na}=3.0\times 10^{22}\mathrm{cm}^3`$ denote the mass and polarizability of the sodium atom, respectively. (Note that polarizability in the Lorentz-Heaviside which we use is $`4\pi `$ times that in Gaussian units.) If $`t`$ is of order $`10^3s`$ (the time needed for an atom with a kinetic energy of order 300K to travel a few centimeters), and $`z`$ is of order $`1\mu \mathrm{m}`$, the fractional deflection is significant.

An alternative to measuring the deflection of the beam might be to measure the relative phase shift in an atom interferometer, which are sensitive to phase shifts of the order of $`10^4`$ radians WPW . If one path of the interferometer were to be parallel to the focal line at a mean distance of $`a`$ and the other path far away, the phase difference will be

$$\mathrm{\Delta }\varphi =\frac{t}{2}\alpha 𝐄^2=0.04\left(\frac{\mathrm{\Lambda }}{10^3}\right)\left(\frac{\alpha }{\alpha _{Na}}\right)\left(\frac{1\mu \mathrm{m}}{a}\right)^4\left(\frac{t}{10^3s}\right).$$
(65)
Both deflection and phase shift measurements would require the atoms to be rather well collimated and localized near the focal line.

Perhaps the most dramatic confirmation of enhanced vacuum fluctuations would be the trapping of atoms near the focus. This would require a mirror with $`\theta _0>\frac{\pi }{2}`$ so that there is a region with $`𝐄^2>0`$. It would also require the atoms to be cooled below a temperature of about

$$T_m=2\times 10^9K\left(\frac{\mathrm{\Lambda }}{10^3}\right)\left(\frac{\alpha }{\alpha _{Na}}\right)\left(\frac{1\mu \mathrm{m}}{a}\right)^4.$$
(66)
Thus at temperatures of the order of $`10^9K`$, atoms could become trapped in a region of the order of $`1\mu \mathrm{m}`$ from the focus. This type of trapping would be quite different from that currently employed Phillips  in that it would require no applied classical electromagnetic fields.

The peculiar property that the force on an atom can be attractive from certain directions and repulsive from other directions requires some discussion. If one approaches from an attractive direction, the potential energy is becoming increasingly negative, $`V1/a^4`$, whereas from a repulsive direction it is becoming large and positive, $`V+1/a^4`$. In both cases, $`V0`$ for large $`a`$. At some minimum value of $`a`$, the inverse fourth power behavior has to be modified. Thus the exact potential energy $`V(a,\gamma )`$ must be a continuous function with a stable minimum, as illustrated in Fig. 13. Note that the requirement that $`V0`$ for large $`a`$ in all directions rules out the possibility that there is only a saddle point, and no minimum. Although we can infer the existence of the minimum from our results, the detailed form of $`V(a,\gamma )`$ near the minimum could only be found by an analysis which goes beyond the approximations made in the present paper.

## IX Discussion and Conclusions

In this paper we have further developed a geometric optics approach to the study of vacuum fluctuations near the focus of parabolic mirrors. The main result of this approach is that the mean squared scalar field $`\phi ^2`$ and mean squared electric field $`𝐄^2`$ grow as inverse powers of $`a`$, the distance from the focus, for small $`a`$. The key justification of geometric optics is its self-consistency. When $`\phi ^2`$ and $`𝐄^2`$ are large, the dominant contributions must come from short wavelengths for which geometric optics is a good approximation. This was discussed more quantitatively in Sect. VII.

We have given some explicit analytic and numerical results for $`\phi ^2`$ and $`𝐄^2`$ near the focal line of a parabolic cylindrical mirror. In this paper, we restricted our attention to the case that the angular size of the mirror, $`2\theta _0`$, is less than $`\frac{4\pi }{3}`$. This insures that there are never more than two reflected rays for a given incident ray and simplifies the analysis We find that for smaller mirrors, $`\theta _0<\frac{\pi }{2}`$, $`\phi ^2>0`$ and $`𝐄^2<0`$ everywhere. In this case, an atom will feel a repulsive force away from the focus. For larger mirrors, $`\frac{\pi }{2}<\theta _0<\frac{2\pi }{3}`$, these mean squared quantities can have either sign, depending upon the direction from the focus. In directions nearly perpendicular to the symmetry axis, $`\gamma \frac{\pi }{2}`$, we find $`\phi ^2<0`$ and $`𝐄^2>0`$. In this case, the force on an atom is attractive toward the focus and trapping becomes a possibility.

In the geometric optics approximation, the mean squared electric and magnetic fields are equal, so the local energy density is equal to $`𝐄^2`$. Thus, when $`𝐄^2<0`$, the local energy density is negative, and one has columns of negative energy density running parallel to the focal line of the parabolic cylinder.

The feasibility of experiments to observe the effects on atoms near the focus was discussed in Sect. VIII. Although the effects are small, it seems plausible that they could be observed.

###### Acknowledgements.

We would like to thank Ken Olum for valuable discussions. This work was supported in part by the National Science Foundation under Grant PHY-9800965, by Conselho Nacional de Desenvolvimento Cientifico e Tecnologico do Brasil (CNPq), and by the U.S.Department of Energy (D.O.E.) under cooperative research agreement DF-FC02-94ER40810.

*

## Appendix A

In this Appendix, we will develop a method which goes beyond the geometric optics approximation, and use it to estimate the size of the corrections due to diffraction effects. Here we will discuss only a scalar field which satisfies Dirichlet boundary conditions on the mirror. The basic method is to write down and then approximately solve an integral equation for the scattered wave. This method was first used by Kirchhoff to discuss diffraction in classical optics. Let $`\varphi (𝐱)`$ be a solution of the Helmholtz equation

$$^2\varphi +k^2\varphi =0,$$
(67)
and let $`G(𝐱,𝐱^{})`$ be a Green’s function for this equation

$$^2G+k^2G=\delta (𝐱𝐱^{}).$$
(68)
If we multiply Eq. (67) by $`G`$, Eq. (68) by $`\varphi `$, take the difference, and integrate over an arbitrary spatial volume $`V`$, the result is

$$\varphi (𝐱^{})=_S𝑑𝐚[\varphi (𝐱)G(𝐱,𝐱^{})G(𝐱,𝐱^{})\varphi (𝐱)].$$
(69)
Here $`S`$ is the boundary of $`V`$, and $`d𝐚`$ is outward directed.

We are interested in calculating the scattered wave from a open surface $`S_1`$, rather than the solution on the interior of a closed surface $`S`$. However, we can let the closed surface $`S`$ consist of three segments, as illustrated in Fig. 14. The first is the surface of interest $`S_1`$, the second is a segment $`S_2`$ which is hidden from the view of an observer in the region where we wish to find the scattered wave, and the third $`S_3`$ closes the surface at a large distance. We now ignore the contributions of $`S_2`$ and $`S_3`$, and assume that $`\varphi =0`$ on $`S_1`$, so that we have

$$\varphi (𝐱^{})=_{S_1}𝑑𝐚\varphi (𝐱)G(𝐱,𝐱^{}).$$
(70)
This is still an integral equation which relates values of $`\varphi `$ on the surface to those off of the surface. We can solve it in an approximation in which an incident wave scatters only once from the surface. Let $`\varphi _{in}`$ be the incident wave, and assume that on $`S_1`$

$$\widehat{𝐧}\varphi =2\widehat{𝐧}\varphi _{in}.$$
(71)
We next choose the empty space Green’s function

$$G(𝐱,𝐱^{})=\frac{\mathrm{e}^{ik|𝐱𝐱^{}|}}{4\pi |𝐱𝐱^{}|},$$
(72)
and interpret Eq. (70) as giving the scattered wave, $`\varphi _s`$,

$$\varphi _s(𝐱^{})=2_{S_1}𝑑𝐚\varphi _{in}(𝐱)G(𝐱,𝐱^{}).$$
(73)
The physical interpretation of this expression is that each point on $`S_1`$ radiates a scattered wave proportional to $`\widehat{𝐧}\varphi _{in}`$; the superposition of these individual contributions forms the net scattered wave, as expected from Huygen’s principle.

Let the incident wave be a plane wave

$$\varphi _{in}=\mathrm{e}^{i𝐤𝐱},$$
(74)
and the surface $`S_1`$ be the portion of the $`x=0`$ plane in the interval $`y_0<y<y_0`$. That is, it is a strip which is infinite in the $`z`$-direction. Take the wavevector of the incident wave to be $`𝐤=k(\mathrm{cos}\theta ,\mathrm{sin}\theta ,0)`$ and the observation point to be $`𝐱^{}=(b,0,0)`$, that is, at a distance $`b`$ from the mirror, as illustrated in Fig. 15. We can now write the scattered wave as

$$\varphi _s(𝐱^{})=\frac{ik\mathrm{cos}\theta }{2\pi }_{y_0}^{y_0}𝑑y\mathrm{e}^{ik\mathrm{sin}\theta y}_{\mathrm{}}^{\mathrm{}}𝑑z\frac{\mathrm{e}^{ik\sqrt{y^2+z^2+b^2}}}{\sqrt{y^2+z^2+b^2}}.$$
(75)
The $`z`$-integration may be performed explicitly, with the result

$$\varphi _s(𝐱^{})=\frac{1}{2}k\mathrm{cos}\theta _{y_0}^{y_0}𝑑y\mathrm{e}^{ik\mathrm{sin}\theta y}H_0^{(1)}(k\sqrt{y^2+b^2}),$$
(76)
where $`H_0^{(1)}`$ is the Hankel function of the first kind. In general, the remaining integral cannot be evaluated explicitly. However, in the case of an infinite plane mirror, $`y_0\mathrm{}`$, it can be evaluated, with the result

$$\varphi _s(𝐱^{})=\mathrm{e}^{ik\mathrm{sin}\theta b}=\mathrm{e}^{i𝐤^{}𝐱^{}}.$$
(77)
This is just the result predicted by geometric optics. In the special case of an infinite plane mirror, geometric optics gives the exact result. In the case of a finite mirror, Eq. (77) is still a good approximation in the high frequency limit. In this limit, one may evaluate the integral in Eq. (76) using the stationary phase approximation. So long as the classical path intersects the mirror in the interval $`y_0<y<y_0`$, then there is one point of stationary phase within the range of integration. This intersection occurs if $`y_0<b\mathrm{cos}\theta <y_0`$. The result of the stationary phase approximation is just Eq. (77), reflecting the fact that geometric optics is a good approximation in the high frequency limit.

Now we wish to give a quantitative estimate of the accuracy of the approximation. Let

$$\mathrm{\Delta }\varphi _s=\frac{1}{2}k\mathrm{cos}\theta \left(_{\mathrm{}}^{y_0}+_{y_0}^{\mathrm{}}\right)dy\mathrm{e}^{ik\mathrm{sin}\theta y}H_0^{(1)}(k\sqrt{y^2+b^2}).$$
(78)
This is the difference between the stationary phase (geometric optics) approximation and the exact result for a finite mirror, so $`|\mathrm{\Delta }\varphi _s|=|\mathrm{\Delta }\varphi _s|/|\varphi _s|`$ is a fractional measure of the accuracy of the approximation.

In the high frequency limit that $`kb`$, we can use the large argument form for $`H_0^{(1)}`$,

$$H_0^{(1)}(z)\sqrt{\frac{2}{\pi z}}\mathrm{e}^{i(z\frac{\pi }{4})},$$
(79)
and write

$$\mathrm{\Delta }\varphi _s\sqrt{\frac{k}{2\pi }}\mathrm{cos}\theta \mathrm{e}^{i\pi /4}_{y_0}^{\mathrm{}}𝑑y(y^2+b^2)^{\frac{1}{4}}\mathrm{e}^{ik(\mathrm{sin}\theta y\sqrt{y^2+b^2})}+(y_0y_0).$$
(80)
The above integral still cannot be evaluated explicitly, but we can estimate it as being of order $`1/(k\sqrt{y_0})`$ when $`ky_0`$. With this estimate, we find that

$$|\mathrm{\Delta }\varphi _s|O\left(\frac{\lambda }{y_0}\right)^{\frac{1}{2}},$$
(81)
where $`\lambda `$ is the wavelength of the incident wave.