1. Introduction

Conical refraction is arguably one of the earliest theoretical predictions of science. In 1832, Hamilton[1], whose works have had a fundamental impact on mechanics and optics and their interrelation, predicted that an incident light beam traveling along one of the optic axes of a biaxial crystal would produce a whole cone of light rather than separate ordinary and extraordinary refracted rays. The prediction was experimentally verified shortly thereafter by Lloyd [2]. The reason is the nature of the dispersion relation F(ω,k_x ,k_y ,k_z ) = 0. The constant frequency surfaces are quartic surfaces in wave vector (k_x ,k_y ,k_z ) space which have the special property that the two sheets meet at four points. The lines joining opposite pairs are the optic axes. This unusual feature gave rise to studies of a general class of surfaces known as Kummer surfaces whose branches meet at points. At these points (called diabolical points by Michael Berry), the two touching surfaces are cones. For small differences in the refractive indices in the three directions x, y, z (for the photo-refractive crystals KTP, LBO, Sn₂P₂S₆ and KNbO₃ respectively at a wavelength of 0.5μm, the ratios (n_x : n_y : n_z ) are (1.78:1.79:1.89), (1.58:1.61:1.65), (3.05:3.17:3.29) and (2.21:2.33:2.40), the cone angle is very wide. The dual cone of the pencil of normals along which the light rays travel has a narrow angle which we call 2p. Therefore, each ray in the cone is almost parallel to the optic axis. Because of this one can treat the complex amplitudes A and B of the two electric field polarizations as envelopes which vary slowly in space and time and the electric field itself as a product of these amplitudes and a carrier wave propagating along the optic axis. The behavior of a conically refracted beam in a nonlinear medium with diffraction can therefore conveniently be described by the solutions of two coupled partial differential equations of nonlinear Schrödinger type. The equations, which we call the CR equations, are new, and their solutions exhibit some surprising and unexpected behaviors.

Our results are consistent with the analysis and experiments of Berry [3], and Berry, Jeffrey and Lunney [4] who have analyzed the influence of diffraction on the Poggendorff outer and inner bright and intermediate dark rings and far field central spike, a manifestation of the fact that the Fourier transform of the incident beam has a large zero wavenumber component. Berry’s analysis allows for a finite length (L) crystal (after which the cone becomes a cylinder)) and captures in exquisite detail the diffraction induced decorations on the rings and central spike in the limit where the ratio of pL/w ₀ (w ₀ is the beam width) is large, namely when the radii of the rings is much larger than the beam width. The present letter explores the competing influences of conicity, diffraction and nonlinearity. We announce two fundamentally new and qualitatively different features. First, in the nonlinear system, plane waves are modulationally unstable even in the case where the nonlinearity is defocussing. The connection with the modulational instability of coupled beams in defocussing media [5, 6, 7] is discussed. Second, there is an interesting analogue to Helmholtz-Gauss non-diffracting beams [8, 9, 10] in that, in the absence of diffraction and nonlinearity, certain incident beam profiles undergo no conical refraction. If the incident beam is analytic, namely if the the combination B + iA = F(ζ) is analytic in ζ = x + iy (x and y are the coordinates perpendicular to the direction of the optical axis), then the beam displays no conical spreading. The simplest example is an optical vortex F(ζ) = ζ. Moreover, an incident beam which is a modulation of a Gaussian beam G(r) by an analytic F(ζ) (e.g. a Gaussian beam with a vortex singularity) will maintain the vortex near the optic axis. The conical rings emanating from the Gaussian beam part are modified in nontrivial ways.

Our goal in this short letter is to (a) introduce the CR equations and discuss briefly how they are derived and (b) discuss their most interesting properties. A second publication, which will include the results of experiments being carried out by Residori and Bortolozzo at INLN, will provide more details and results.

2. Analysis

From Maxwell’s equations ∇×H= ∂D/∂t, ∇ × E = -∂B/∂t, ∇ · D = ∇ · B = 0, with D = ε ₀ ($n_x^2$ E_x , $n_y^2$ E_y , $n_z^2$ E_z ) +P_NL , and $n_x^2$<$n_y^2$=n ² <$n_z^2$, we first seek linear monochromatic solutions with E⃗ = $\widehat{E}{e}^{i\overrightarrow{k}\cdot \overrightarrow{x}-i\omega t}$ and P⃗_NL = 0 which yields ΩÊ = 0. The real symmetric matrix is

Setting det Ω = 0 and fixing ω, we obtain a quartic surface whose optic axis is given by $\overrightarrow{k}=\frac{\mathit{n\omega }}{c}$(sinθ,0,cosθ) where $\sin \theta =\frac{n_z}{n}\sqrt{\frac{n^2-n_x^2}{n_z^2-n_x^2}},\cos \theta =\frac{n_x}{n}\sqrt{\frac{n_z^2-n^2}{n_z^2-n_x^2}}.$ At this k⃗ value, ω has two (left and right) eigenvectors ($(\frac{1}{n_x^2}\cos \theta ,0,\frac{-1}{n_z^2}\sin \theta )$) and (0,1,0). We now seek solutions of the spatially modulated and nonlinear system in the form

The complex amplitudes B and A are slowly varying functions of x, y, z and t so that $\frac{c}{\mathit{n\omega }}\nabla B\mathrm{and}\frac{1}{\omega }\frac{\partial B}{\partial t}$ are much smaller than B. The correction E⃗ ₁ is also small. Small here is given by a measure of the semi-angle of the cone normal to the constant frequency surface at the diabolical point $p=\frac{1}{n}\sqrt{\left(n-n_x\right)\left(n_z-n\right)},$ typically between 0.03 and 0.1. Note that to leading order, ∇·D⃗ = ∇ · ($n_x^2$ E_x , n ² E_y , $n_z^2$ E_z ) = 0. Indeed, it is more natural to deal with the displacement vector D⃗ (its linear approximation) because its two polarization directions and k⃗ form a right angled triad coordinate system, z′ = xsinθ+zcos θ,x′ =xcosθ - zsinθ andy′ =y.

The long distance and time evolution of B and A are obtained in the standard way. We insert E⃗ as given in equation (2) into the Maxwell equations. The slow derivatives are incorporated by setting $\mathit{\omega \to \omega }+i\frac{\partial }{\partial t},\overrightarrow{k}\to \overrightarrow{k}-i\nabla$in Ω and expanding in a Taylor series treating $\frac{1}{\omega }\frac{\partial }{\partial \mathit{t}}\mathrm{and}\frac{c}{\mathit{n\omega }}\nabla$ as small. The nonlinear terms, arising from P⃗_NL are easily added. The equation for the correction E⃗ ₁ in (2) has on its right hand side vectors proportional to exp $\exp \mathit{\pm i}\left(\frac{\mathit{n\omega }}{c}\left(x\sin \theta +z\cos \theta \right)-\mathit{\omega t}\right)$ and higher harmonics. Solvability demands that the vector multiplying the first harmonic has no projection onto the eigenvectors of Ω. This gives two equations for B and A which we call the equations for nonlinear conical refraction, the CR equations. We write B and A as functions of x′ = xcosθ - zsinθ, y′ = y, z′ = xsinθ +zcosθ and drop the primes. After considerable manipulation, we obtain the beautiful and simple equations

The parameter $\beta =\frac{c}{2\mathit{n\omega }},$ and the Laplacian ∇² is the two dimensional one. Its simplicity depends on the fact that p is small. This also means that the three terms, conicity, diffraction and nonlinearity can balance. For wide cone angles, the diffraction terms are much more complicated and conicity dominates. The different values of the nonlinear coefficients τ ₁ and τ ₂ and the non-unit cross coupling coefficient γ reflects the fact that many crystals are also nonlinearly non-isotropic. In orthorhombic crystals, both electro-optic and photovoltaic effects are important. The different electro-optic coefficients determine the degree of anisotropy, the magnitudes of the nonlinear coefficients τ ₁, and τ ₂ and the cross coupling coefficient γ. The nonlinear refractive index (|A|² + γ|B|²)(1 + δ (|A|² + |B|²))^-1 is due to the photovoltaic field induced along the z-axis. For KNbO₃, τ ₂ may be as much as 6τ ₁, and γ may be as small as 0.2. The saturation intensity δ ^-1 is determined by the dark intensity [10]. For the purposes of this letter, we consider the isotropic, Kerr nonlinearity case where τ ₁ = τ ₂, γ= 1 and δ = 0. Different values will be treated later. We will also deal with continuous monochromatic incident beams so that $\frac{\partial A}{\partial t}=\frac{\partial B}{\partial t}=0.$ Pulsed beams, however, can easily be handled.

3. Results

A comprehensive analytical, numerical and experimental study of conical refraction with non-linearity will appear later. In this letter format, we present, in list form, some of the main results.

1. Poggendorff rings and the Berry analysis.

Solving equation (3) with β = τ = 0 and B(x,y,0) = iA(x,y,0) = e ^{-r2/2$w_0^2$}, where r ₂ = x ² + y ², we find that B(x,y, z) and A(x,y, z) can be expressed as Fourier integrals with exponential phase h(ℓ,m) = -$w_0^2$ k ²/2ℓipkz + ilx + imy, k⃗ = (l,m), $k=\sqrt{{\ell }^2+m^2}.$ Standard steepest descent methods show that the integrals are dominated by the choices ℓ/k = x/r, m/k =y/r (consistent with the ray theory) and k = ±i(r∓pz)/$w_0^2$. At these points the absolute value of the determinant of the matrix of second derivatives of h is given by $w_0^2$ r/k whose inverse is a measure of the weight of the contribution from the given k = ±i(r∓pz)/$w_0^2$. This leads to intensities proportional to |r ∓pz|e ^{-(r±pz)2/$w_0^2$}. The minus sign contribution dominates. The intensity I has two bright peaks at either side of the dark ring. These are the Poggendorff rings and the zero arises because the weight of the light emanating from the conical point r = 0 at z = 0 is zero. Berry’s results can be obtained by solving equation (3) with τ = 0 and finite β. Indeed, the Fourier integrals are exactly those obtained in Berry[3].

2. No y dependence, no conical refraction!

For beams with no dependence on the y, the direction of polarization of the ordinary wave B there is only nonlinear cross coupling of the ordinary and the extraordinary beams. In this case, there is no conical refraction. However (see note 5 below), these solutions are nonlinearly unstable to y dependent perturbations.

3. Initially CR (Cauchy-Riemann) means no CR (Conical refraction)!

For β = τ = 0, we recognize the right hand side of the CR equations (3) as the Cauchy-Riemann conditions (consistent with the fact that both B and A satisfy a wave equation $\left(\frac{{\partial }^2}{{\partial z}^2}-p^2{\nabla }^2\right)\left(B,A\right)=0).$ If we let the subscript j denote taking either the real or the imaginary component of a complex amplitude, then if B_j + iA_j is analytic in ζ = x + iy at z = 0, there is no conical refraction.. We will discuss below what happens when at z = 0 the electric field is circularly polarized with B = iA by writing B+iA = F(ζ)F ₊ (r, α,z), and B - iA=F(ζ)e^iα F _- (r, α,z) where F(ζ) = ζ is an optical vortex and F _±(r, α,0) are Gaussian beams. In a later paper, we will also discuss the effects of nonlinearity and conical refraction in incident diffractionless beams [8, 9] and explore whether there are beams which neither diffract nor conically refract.

4. Linear dispersive relation.

Linear waves with A, B proportional to exp(iℓx + imy - ikz) have the dispersion relations k _± = βk ² ±pk, k ² = ℓ² + m ², with different combinations of A and B taking on the two signs. For long waves, pk ≫βk ², conical refraction waves propagate like sound waves. For χ ⁽²⁾ non-linearities, it would be interesting to explore their wave turbulence properties. The reason is as follows. Pure sound waves share energy via triad interactions which redistribute energy along rays in wavevector space. It is not known how kinetic energy is then redistributed between different rays [11]. The CR equations with χ ⁽²⁾ nonlinearities would make an excellent vehicle for examining this question.

5. The modulational instability.

Perturbing a right circularly polarized plane wave B(z) = iA(z) = B ₀ e ^iτI₀z, I ₀ = 2|B ₀|² in the form ΔB = (be ^{σz+lℓx+lmy} e ^iτI₀z, ΔA = (ae ^{σz+iℓx+imy}) e ^iτI₀z, we find

 

Fig. 1. Evolution of amplitude of the most unstable mode, semi-log scale. The x’s are for a perturbation with k = 11 and τI ₀ = 2, the o’s are for k = 9 and τI ₀ = -2. The lines are for reference, and show the predicted growth rates. Parameter are β = 0.1, p = 1

Download Full Size | PDF

where k ² = ℓ² + m ². If p = 0, we obtain the usual Benjamin-Feir result that instability occurs in the window 0 < β ² k ² < 2βτI ₀ and thus only if the product βτ is positive. Since β = c/2nω > 0, instability only occurs if the nonlinear refractive index τ is positive, the focussing case. For nonzero p, there is always instability. For βτ > 0, it occurs in the short wave window p ² < β ² k ² < p ² + 2βτI ₀. For βτ < 0, it occurs in the window min(p ² + 2βτI ₀, 0) < β ² k ² < p ². Therefore, plane waves are always unstable. The combination of nonlinearity and conicity always produces a window of instability in the wavenumber. The window is absent if there is no nonlinearity. Nonlinearity makes conical refraction inevitable. For p ² ≫ 2βτI ₀, the most unstable wave number is k = p/β + τI_o /2p and the corresponding growth rate is close to τI ₀/2. For other polarization choices, the result is similar. In Fig. 1, we plot the z growth of the amplitude of the most unstable perturbation on a semi-log plot. In these cases the predicted slope is ≈ 1.0945 for τI ₀ = 2 and k=11, and ≈ 0.8955 for τI ₀ = -2 and k = 9 which agrees well with the measured slopes.

One might think that this analysis suggests that in a nonlinear medium the spreading conical rings might destabilize to circumferential perturbations. Two factors mitigate against this. First, the spreading of the rings means that the most unstable wavenumber at z = z ₁, say, slides out of the instability window as z evolves. Second, the amplitudes in the rings decrease and are much smaller than those at the center of the original beam. Numerical simulations verified this conclusion. No evidence of the ring breaking into spots, as would happen for a fixed radius beam in a focussing medium, were found.

This result is technically similar to the modulational instability found by Roskes [5] and Wright [6] for coupled beams in defocussing media. In that case, the key feature is that the two beams have slightly different wavenumbers so that the beating of one beam against the other provides for the destabilization of the plane wave. Here the detuning effect is captured by p. To see this, take (3) with p = 0 and $p=0\mathrm{and}\frac{\partial }{\partial y}=0$ (as is done in [6]) and then look for plane wave solutions with slightly different wavenumbers (B,A) = e ^{±iqx/2β-iq2z4β}(b,a). The variables b and a satisfy (3) with q = p and $q=p\mathrm{and}\frac{\partial }{\partial y}=0.$ For CR propagation, no detuning is required. It is automatically provided by conical refraction. It is important, however, to recognize that the plane waves are also unstable to purely y dependent perturbations.

6. Typical parameter values.

If w ₀ is the incident beam waist, the conical, diffractive and nonlinear lengths ℓ_c , ℓ_d and ℓ_n are respectively ${\ell }_c=\frac{w_0}{p},{\ell }_d=\frac{w_0^2}{\beta }\mathrm{and}{\ell }_n={\left(\tau I_0\right)}^{-1}.$ For currently available crystals, the values are respectively (p = 0.04, for λ = 0.5μm, β = 2 × 10^-8m, τI ₀ ≈ 2 × 10²m^-1, w ₀≈2× 10^-3m) ℓ_c ≈ 50mm, ℓ_d ≈ 100m, ℓ_n ≈ 10^-2m although the diffractive distances are shortened considerably by the local field gradients which can be much larger than $w_0^{-1}$ . Unfortunately, from the experimental viewpoint, pure biaxial nonlinear crystals such as KNbO₃, which are currently available are small, typically less than 1cm in diameter. We will have to await larger crystals or perhaps natural contexts, such as perhaps the interstellar medium which in certain situations might behave as a biaxial crystal and/or larger nonlinearities before some of the predictions of this letter are experimentally verifiable. We can, however, verify the predictions using numerical simulations.

 

Fig. 2. Total intensity at z = 1 for (a) Gaussian with τ = 0, (b) with τ = 1, (c) modulated Gaussian (F(ζ) = ζ) with τ = 0, and (d) with τ = 1.

Download Full Size | PDF

7. Gaussian beams with and without analytic modulation.

It is convenient to rewrite equation (3) in polar coordinates. The angular dependence of the coefficients can be removed by (B,A) ^T = C(f,g) ^T where C is the matrix which rotates vectors clockwise by α/2. This transformation reveals that the electric field polarization reverses on completing one 2π cycle. Since E⃗ and A,B are single valued, this twist is compensated by an opposite rotation of the phases of A and B. For example, for an initially circularly polarized beam with A = -iB, f and g are proportional to e ^-iα/2.

8. A convenient transformation to study analytic modulation of Gaussian beams

For initially circularly polarized beams, make the transformation B + iA = F(ζ)F ₊(r,α,z), B - iA = e ^iα F(ζ)F _-(r,α,z) where ζ = x + iy and F(ζ) is analytic in ζ. Then equation (3) becomes

For $F\left(\zeta \right)=c{\zeta }^s,\frac{2F\prime e^{\mathit{i\alpha }}}{F}=\frac{2s}{r},\mathrm{and}\mathit{FF}*={\mid c\mid }^2{r}^{2s}.$ For initial conditions for which B + iA = F(ζ) and B - iA = 0, F _-, representing the opposite circular polarization, remains zero, and F ₊ develops a phase change. But until the modulational instability develops, there is no conical refraction.

9. Numerical simulation of modulated beams

We present in Figs. 2 and 3 comparisons of the propagated beam in the cases of (a) a Gaussian beam with no nonlinearity (τ = 0) and no optical vortex (F(ζ) = 1), (b) a Gaussian beam with nonlinearity (τ = 1) and no optical vortex (F(ζ) = 1), (c) a beam with no nonlinearity (τ = 0) and a modulating optical vortex (F(ζ) = ζ), (d) a beam with nonlinearity (τ = 1) and a modulating optical vortex (F(ζ) = ζ). The parameter values used were p = 1, β = 0.1 and w = 1. In Fig. 2 the profiles are shown for beams at z = 1. In Fig. 3 the beam intensities are plotted as a function of r = √x ² + y ² at z = 0 (black), z = 1 (blue) and at z = 2 (red). The computational domain has width 6π. We observe the following. In Figs. 2a and 3a the rings are just beginning to form. They are not clearly visible in Fig. 2a (at z = 1), but the intensity profiles in Fig. 3a have acquired local maxima. Note that at z = 1, the central spike has not yet formed, but it has begun to form at z = 2 as the red curve shows. Adding nonlinearity (as in Figs. 2b and 3b) dramatically increases the central intensity and the intensities of the rings. The rings are sharper and thinner and have smaller radius due to the nonlinear focussing effect. In Figs. 2c and 3c, an optical vortex has been added. It maintains its identity as it propagates while the Gaussian part of the beam spreads. In Figs. 2d and 3d, nonlinearity has been added. This time, however, the intensity deficit at the vortex center allows the rings to spread further than in case (b).

 

Fig. 3. Total intensity graphed radially at z=0 (black), 1 (blue), and 2 (red) cases as in Fig. 2 The peak for (b) at z = 2, r ≈ .2 goes off scale to ≈ 6.8

Download Full Size | PDF