Nonlinear optics utilizing composite optical material architectures for providing enhanced nonlinear susceptibility
DESCRIPTION 
     The present invention relates to the field of nonlinear optics and particularly to optical devices, systems and methods of processing optical radiation which rely upon nonlinear optical properties of a body of composite material having two or more components which responds to an applied optical field due to incident optical radiation nonlinearly with the strength of the optical field. 
     The invention is especially suitable for providing optical systems, devices and optical signal processing methods which rely on the nonlinear response of a material system. The nonlinear response may affect polarization of an optical beam, the amplitude of the beam (in general an optical signal) which is incident upon the material system and which propagates therein. The nonlinear response may change the index of refraction of the material system and produce other nonlinear optical interactions. Interactions and processes, as may be provided by methods and systems embodying the invention include harmonic generation, such as doubling and tripling of the frequency of the incident radiation, mixing of two beams for sum frequency generation and difference frequency generation, optical parametric oscillation, nonlinear or saturable absorption and photonic (including electro-optic effect dependent processes). The term optics or optical as used herein includes light radiation in the visible, the infrared and ultraviolet region. Aspects of the invention may be used in effecting nonlinear processes above the ultraviolet and below the infrared range, for example x-ray radiation. The invention is believed to be primarily of interest in connection with optical radiation in the visible, near infrared and ultraviolet frequency ranges. 
     The present invention provides improvements in nonlinear optical processing through the use of a body of material (a material system) having an architecture where different materials are combined to form a composite. Preferably the composite is in the solid state, but it may be at least partially in the liquid state and may even be partially in the gaseous state, when the components in the gas are physically stable. There are at least two components, one of which is the host and the other of which provides inclusions in the host. The host is of one material type and the components may differ in their type of material, and significantly in their refractive indices from the index of the host. The inclusions are particles which are dispersed in the host so that the material system is homogeneous on a distance scale of the order of the wavelength of the radiation being processed in the material system. The host and inclusion materials are described by means of their linear or nonlinear responses which depend on the intensity of the radiation propagating in the material system. These responses are more precisely defined by the effect of linear and nonlinear optical susceptibilities of the materials. 
     It is the principal feature of this invention to provide improved systems, devices and methods for processing of optical radiation, which embody composite material systems presenting nonlinear susceptibility larger than of any component of the composite above. Such increased or enhanced nonlinear susceptibility increases the efficiency of nonlinear optical interactions, thereby improving the performance of systems and processes depending upon such interactions. Moreover, the invention permits enhancement of nonlinear optical properties through the use of materials, both for the host and inclusions in the composite system, which are optically transmissive and do not absorb the radiation which propagates therein, and since the radiation intensity is not depleted by absorption, the processes are carried out more efficiently than where absorbing materials are needed. 
     Composite material systems have been studied experimentally and theoretically. The pioneer work was that of J. C. Maxwell Garnett, which appeared in Philos. Trans. Roy. Soc., London 203, 385 (1904) and 205, 237 (1906). The Garnett work considered only liquid or glass hosts with metallic particles. Both the host and the particles were characterized as having linear responses. The effects obtained (color variations) depended upon absorption. Nonlinear optics was not involved. The field of nonlinear optics did not exist until the development of the laser which provided radiation of sufficient intensity to excite nonlinear effects. Composite systems involving nonlinear particles in a linear host have been reported. See D. Ricard, P. Roussignol and C. Flytzanis, Optics Lett. 10, 511 (1985); F. Hache, D. Ricard and C. Flytzanis, J. Optic Soc. Am., B 3,1647 (1986); K. C. Rustagi and C. Flytzanis, Opt, Lett. 9, 334 (1984); P. Roussignol, D. Ricard, J. Lukasik and C. Flytzanis, J. Opt. Soc. Am. B 4, 5 (1987); G. S. Argarwal and S. Dutta Gupta, Phys. Rev. A 38,5678 (1988); and J. W. Haus et al., J. Opt. Soc. Am. B 6,797 (1989). All of the material systems reported utilized a host (water or glass) with a linear, rather than a nonlinear response. Nonlinear, but absorbing, particles were reported in the articles authorized by Flytzanis and others and enhancements in the value of the nonlinear susceptibility were reported. Such enhancements in nonlinear response were measured against either the nonlinear susceptibility of the components of the composite or a nonlinear material which has often been used as a standard, namely carbon disulfide (CS 2 ). 
     The present invention is based upon the discovery that the presence of a nonlinear host in the composite (more than a linear host which merely supports or suspends the nonlinear particles) makes an appreciable contribution to the nonlinear susceptibility of the composite. Moreover, the use of a nonlinear host enables the enhancement of the optical nonlinearity to be achieved even where the materials are nonabsorbing (for the case in which both components are lossless or at least the host is lossless). By sacrificing, some absorption in the inclusion component, even greater enhancement in optical nonlinearity can be achieved. 
     As noted above, the earliest theory to deal with the optical properties of such composites is due to Maxwell Garnett. There the composite material is assumed to be comprised of spherical inclusion particles embedded in a host material. Both the inclusion particles and host are isotropic and respond linearly to the incident light. The article by Argarwal and Dutta Gupta, and the article by Haus et al., referred to above, present theoretical studies of composite nonlinear optical materials based on the Maxwell Garnett model with a nonlinear response in the inclusion material, but not in the host. 
     We present a full generalization of the Maxwell Garnett model which shows the nonlinear correction to the effective medium dielectric constant and consider the case where either or both the host and inclusion components possess a third-order nonlinear optical response, which for isotropic materials can be characterized by two parameters usually denoted A and B, and in terms of those parameters and the dielectric constants of the constituent materials we derive expressions for the parameters A and B of the effective medium. The parameters A and B express the third-order nonlinear susceptibility as explained in Section 4.2 of the text, Nonlinear Optics, by Robert W. Boyd, Academic Press (1992). We show that the composite materials can possess a nonlinear susceptibility larger than that of its constituents, and for a nonlinear host material, the tensor properties of the nonlinear susceptibility of the composite are considerably enhanced (by a factor of five) over the pure host material. A similar enhancement, perhaps somewhat smaller, exists for the second-order nonlinear susceptibility. The host is then anisotropic, Material systems which presently are believed useful are also discussed. 
     We present our showing with the following organization. First we present the model of tlie composite material topology and the basic equations that govern the electromagnetic fields in it. Then we give a heuristic derivation of the relation between the mesoscopic and macroscopic electric fields. A more rigorous derivation of essentially the same result is presented in Appendix I. Then we derive the Maxwell Garnett result for the effective-medium dielectric constant. We then consider the case of nonlinearity first in the inclusion material, and then in the host material. Some geometrical formulas used in the section on the case where the host is nonlinear are derived in Appendix II. Example results for composite nonlinear susceptibilities are given in the sections on the nonlinear inclusions and host. We summarize our results and consider the case of nonlinearity in both the host and inclusion materials in a final section. 
     While the inclusions are treated in a topology where the inclusions are spherical and distributed on the host, the particles can be elongated (e.g. elliptical). The topology of a composite material may consist of two or more interdispersed components, i.e. the particles being part of a matrix containing the host. 
     It is therefore the principal object of this invention to provide an improved system for processing of optical radiation utilizing a composite material system presenting an enhanced nonlinear susceptibility to the radiation. 
     It is another object of the invention to provide improved optical devices wherein nonlinear optical interactions occur. 
     It is a still further object of the invention to provide improved methods of processing optical radiation and devices and optical systems, all utilizing composite material systems having enhanced nonlinear susceptibility. 
     It is a still further object of the invention to provide improved systems for optical radiation and signal processing, improved devices in which nonlinear optical interactions occur and improved methods of processing optical radiation and signals nonlinearly, dependent upon the intensity of the radiation or the signals, in which composite material systems having enhanced nonlinear susceptibility are used. 
     Briefly described, therefore, the invention provides a system for processing optical radiation to provide nonlinear response which is a function of the intensity of the radiation. A body of composite material (a material system) is used having a host material component and an inclusion material component. At least the host has an optical response which varies nonlinearly with the amplitude of an optical field applied to the body. Means are provided for applying the field to the body at sufficient intensity to induce the nonlinear response of the body to the field. The nonlinear response is defined by a nonlinear susceptibility greater than the susceptibility of either of the host or inclusion components above. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The foregoing and other objects, features an advantages of the invention as well as the best mode now known for practicing the invention and preferred embodiments thereof will become more apparent from a reading of the following description in connection with the accompanying drawings in which: 
     FIG. 1 is a schematic diagram illustrating the topology of a composite material system (the architecture of a body of composite material) which is provided for use in systems, devices and methods embodying the invention; 
     FIG. 2 is a schematic diagram similar to FIG. 1 showing a sphere of radius R centered on a point r arrow, which may lie either in an inclusion (a particle) or in the host component of the material system; 
     FIG. 3 is a plot showing the variations of the nonlinear coefficients A and B of a composite optical material system with fill fraction f (the fraction of the volume of the material system constituted of the inclusions) for several different values of the ratios of linear dielectric constants; 
     FIG. 4 is a plot showing the enhancement in A+(1/2)B versus the fill factor f of the inclusion material for several values of the ratio ε i  /ε h  (the linear dielective constant of the inclusions ε i  to the linear dielective constant of the host component ε h ) of linear susceptibilities for the case in which the host responds nonlinearly and the inclusions respond linearly; in FIG. 3 and FIG. 4 both the host and material components being nonabsorbing of the optical radiation; and 
     FIGS. 5a, b and c are three plots of the ratio B/A of nonlinear coefficients of the composite material system versus fill factor f of the inclusion component for several values of the ratio ε i  /ε h  of linear dielective constants for the case in which the inclusions respond linearly and the host responds nonlinearly with different host material nonlinear coefficient relationships, B h  =0 in plot (a), B h  =A h  in plot (b) and B h  =6A h  in plot (c). 
    
    
     THE MODEL AND BASIC EQUATIONS 
     We assume a composite topology where small particles, which we call inclusions, are distributed in a host medium. As an approximation, the inclusions are assumed to be spheres of radius a. 
     Preferrably the particles are small with diameters of the order of 50 Angstroms. This reduces scattering since the scattering cross section σ of a small dielectric particle depends on the particle radius a as a 6 , and consequently the attenuaion coefficient α=Nσ, where N=3f/4πa 3  is the number density of inclusion particles, can be made arbitrarily small even at a fill fraction approaching unity by using particles of sufficiently small radius α. For 50-angstrom-diameter inclusion particles, the predicted attenuation coefficeint is less than 1 cm -1 . 
     We define the characteristic distance between inclusion particles to be b, and we assume that 
     
         a&lt;&lt;b&lt;&lt;λ,                                            (2.1) 
    
     where λ is the wavelength of light in vacuum at frequencies of interest (see FIG. 1). Denoting the electric and magnetic field vectors by e(r,t) and b(r,t), respectively, we write the Maxwell equations in the form: ##EQU1## where p(r,t) is the dipole moment per unit volume, and where magnetic effects are assumed to be negligible. We seek stationary solutions of these equations of the form 
     
         e(r,t)=e(r)e.sup.-ωt 30 c.c. .                       (2.3) 
    
     We describe the dielectric response of the material system by means of the relation 
     
         p(r)=χ(r)e(r)+p.sup.NL (r),                            (2.4) 
    
     where ##EQU2## is the spatially varying linear susceptibility, and where p NL  (r) is the non-linear polarization, discussed in detail in the section entitled &#34;nonlinearity in the inclusions.&#34; We define ##EQU3## and combine Eqs. (2.2) through (2.6), the Maxwell equations (2.2) become 
     
         ∇·[ε.sup.h e(r)]=-4π.tbd.·p.sup.s (r) 
    
     
         ∇·b(r)=0 
    
     
         ∇×e(r)-iωb(r)=0 
    
     
         ∇×b(r)+iωε.sup.h e(r)=-4πiωp.sup.s (r),                                                      (2.7) 
    
     with ω≡ω/c, where the inclusion and host linear dielectric constants are given by 
     
         ε.sup.i,h =1+4πχ.sup.i,h,                   (2.8) 
    
     and where we have introduced the &#34;source&#34; polarization 
     
         p.sup.s (r)≡p&#39;(r)+p.sup.NL (r).                      (2.9) 
    
     We have used lower case letters (e, b, etc.) for the electromagnetic fields appearing in these equations because, although we have adopted a macroscopic description of each constituent of the composite, we next introduce a &#34;more macroscopic&#34; description which we obtain by averaging the fields over a volume that contains many inclusions. We perform this average through the use of a smoothly varying weighting function Δ(r)=Δ(r), where r=|r|. The volume integral of Δ(r) is normalized to unity, 
     
         ∫Δ(r-e&#39;)dr&#39;=1,                                  (2.10) 
    
     and it has a range R satisfying the inequalities 
     
         b&gt;&gt;R&gt;&gt;λ.                                            (2.11) 
    
     The averaged fields are then defined by 
     
         E(r)=∫Δ(r-r&#39;)e(r&#39;)dr&#39;, 
    
     
         P(r)=∫Δ(r-r&#39;)p(r&#39;)dr&#39;, 
    
     
         P.sup.NL (r)=∫Δ(r-e&#39;)p.sup.NL (r&#39;)dr&#39;,          (2.12) 
    
     etc., where the integrals are to be performed over all space. It is easily confirmed that this averaging procedure commutes with differentiation, e.g., ##EQU4## where the third of Eqs. (2.13) follows from the second by a partial integration and the assumption that Δ(r-r&#39;)→0 as |r-r&#39;|→∞. By performing an average of the Maxwell equations (2.7), we thus find that 
     
         ∇·[ε.sup.h E(r)]=-4π∇·P.sup.s (r) 
    
     
         ∇·B(r)=0 
    
     
         ∇×E(r)-iωB(r)=0 
    
     
         ∇×B(r)+iωε.sup.h E(r)=-4πiωP.sup.s (r),                                                      (2.14) 
    
     
         where 
    
     
         P.sup.s (r)≡P&#39;(r)+P.sup.NL (r).                      (2.15) 
    
     Our goal is to find a constitutive relation between P s  (r) and E(r). Since the range R of Δ(r) is very much larger than the characteristic separation b of the inclusions. The resulting linear and nonlinear susceptibilities can be expected to be spatially uniform; however, since R&lt;&lt;λ, these averaged fields can be used to describe the propagation of light through the medium (the composite body). The resulting susceptibilities thus characterize an &#34;effective medium&#34; that describes the optical properties of the composite material system, and we refer to the averaged fields E(r), B(r), etc. as &#34;macroscopic fields.&#34; We refer to the fields e(r) and b(r) as &#34;mesoscopic fields&#34;, since they contain more spatial information than the macroscopic fields, but yet are themselves averages of the microscopic electric and magnetic fields, which vary greatly over interatomic distances. 
     MESOSCOPIC AND MACROSCOPIC FIELDS: A SIMPLE PHYSICAL ARGUMENT 
     Since the mesoscopic polarization p(r) is a known function of the mesoscopic electric field e(r) (see Eqs. 2.4 through 2.6 and section 5), we derive a constitutive relation between P(r) and E(r), relating e(r) to E(r). Then pr) can be related to E(r), and the spatial average of p(r) will give P(r) in term of E(r). 
     It turns out that one can derive a relation between e(r) and E(r) that is both rigorous and useful. That relation is derived in Appendix I, and used in the ensuing sections hereof to derive the constitutive relation of the composite medium. In this section we present a simple physical argument that leads to essentially the same result as that derived in Appendix I. This simple physical argument, which follows the original argument of Lorentz, in his text &#34;Theory of Electrons,&#34; Dover, New York (1952). It is somewhat h euristic, but shows why the rigorous result derived in Appendix I is physically reasonable. 
     We wish to relate e(r), which is given by the particular solution of the mesoscopic Maxwell Eqs. (2.7) plus a homogenous solution, to E(r), which is given by the particular solution of the macroscopic Maxwell Eqs. (2.14) plus a homogeneous solution of those equations. We surround the point r by a sphere of radius R centered at r (see FIG. 2), and we write 
     
         e(r)=e.sup.o (r)+e.sup.in (r)+e.sup.out (r) 
    
     
         E(r)=E.sup.o (r)+E.sup.in (r)+E.sup.out (r),               (3.1) 
    
     where e o  (r) and E(r) are the abovementioned homogeneous solutions. We define e in  (r) to be the contribution to e(r) from p s  (r&#39;) taken at points r&#39; within the sphere, and e out  (r) to be the contribution from points r&#39; outside the sphere. Likewise, E in  (r) and E out  (r) contain respectively the contributions to E(r) from P s  (r&#39;) taken at points r&#39; inside and outside the sphere. Now, since E o  (r) is the spatial average of e o  (r) over a distance of the order of R, and e o  (r) varies only over distances of order λ, which is very much greater than R, we may take e o  (r)≃E o  (r). Further, since we have assumed that R&gt;&gt;b, the precise locations of the inclusions outside the sphere are unimportat in determining e out  (r), and to good approximation we may take e out  (r)≃E out  (r). Then, combining Eqs. (3.1), we have 
     
         e(r)-E(r)=e.sup.in (r)-E.sup.in (r).                       (3.2) 
    
     Next, since R&lt;&lt;λ, we can use the laws of electrostatics (the ω→0 limit of the Maxwell equations) in estimating both e in  (r) and E in  (r). For the latter, we ccan also assume that P s  (r is effectively uniform over the sphere, and so Eqs. (2.14) in the electrostatic limit can be shown to yield ##EQU5## We find e in  (r) by solving Eqs. (2.7) in the electrostatic limit, and then restricting the source term p s  (r&#39;) that appears in the solution to points r&#39; within the sphere (|r-r&#39;|&gt;R). The result is can be shown to be ##EQU6## where the radius η of the region excluded from the integration is allowed to go to zero after the integral in Eq. (3.4) is performed. This integral involves the static dipole-dipole coupling tensor for a medium of dielectric constant ε h , which is given by ##EQU7## where r=r/r and U=xx+yy+zz. Defining ##EQU8## and a cut-off function ##EQU9## we can combine Eqs. (3.2), (3.3), and (3.4) to give ##EQU10## This equation is the central result of this section. It shows how the difference between e(r) and E(r) depends only on p s  (r&#39;)at points r&#39; within a distance R&lt;&lt;λof r (cf. Eq. 3.2). This result is used in deriving, to good approximation, a local relation between P s  (r) and E(r) in the ensuing sections hereof. 
     THE MAXWELL GARNETT EQUATIONS 
     In this section we derive the Maxwell Garnett equations, which provide a good approximate description of the linear optical properties of a composite medium of the topology shown in FIG. 1. This result is of course well known (see the articles by J. C. Maxwell Garnett referenced above), but we recover it here both to illustrate our approach and to illustrate the nature of the approximations involved. In addition, we obtain some results that are needed for the calculation of the nonlinear optical response that is presented in the following sections. 
     In both this section and the next, we use the rigorous form of Eq. (3.8), as derived in Appendix I. There, using only the mesoscopic Maxwell Eqs. (2.7), the macroscopic Maxwell Eqs. (2.14), and the definitions (see Eq. (2.12)) of macroscopic fields in terms of mesoscopic fields, we find an exact relation (I.43) between e(r) and E(r) for a Δ(r) characterized by any range R. For R&lt;&lt;λ, that result reduces to ##EQU11## is the so-called &#34;cavity field&#34; (the E(r)-E in  (r) of section 3). Here 
     
         T.sup.c (r)=T.sup.o (r)c(r)                                (4.3) 
    
     is the product of T o  (r) (Eq. 3.6) and a cut-off function c(r) (c(O)=1, c(r)→0 as r∞). The function c(r) is not the simple cut-off function c&#39;(r) that appeared in section 3 (Eq. 3.7), but is given by Eq. (I.60); nonetheless, c(r) has a range on the order of R. 
     We next note three useful results that follow from the properties of the static dipole tensor: ##EQU12## In the first equation, the range of integration is over all space and the result follows from the fact that c(r) depends only on r=|r|, and that the integral of T o  (r) over solid angle, for fixed r, vanishes. In the second two, the range of integration is the interior of a sphere of radius a; in Eq. (4.4b) r is any point within the sphere, while in Eq. (4.4c) r&#39; is a point outside the sphere and r o  the position of the center of the sphere. Eqs. (4.4b) and (4.4c) may be derived, e.g. , by using the fact that T(r) ε h   -1  ∇∇r -1 , and using Gauss&#39; theorem. 
     We now return to a consideration of the polarization of the medium. From Eq. (2.6) we have 
     
         4πp&#39;(r&#34;)=(ε.sup.i -ε.sup.h)θ.sup.i (r&#34;)e(r&#34;),(4.5) 
    
     where we have used Eq. (2.8) and where we define θ i  (r)=1 if the point r&#34; is in an inclusion and θ i  (r&#34;)=0 if r&#34; is in the host. Next, we use Eq. (4.1) and neglect any nonlinear behavior (p NL  (r)=0); then p s  (r)=p&#39;(r) (Eq. 2.9), and Eq. (4.5) yields ##EQU13## where we have put ##EQU14## We now wish to use Eq. (4.6) in the expression for P&#39; see Eq. 2.12), 
     
         P&#39;(r)=∫Δ(r-r)p&#39;(r&#34;)dr&#34;,                         (4.8) 
    
     for which we need p&#39; (r&#34;) at points r&#34; within approximately R of r. We next introduce the standard approximations that are made in deriving the Maxwell Garnett equation. The crucial ansatz is that, at points r&#34; within R of r, p&#39; (r&#34;) within the inclusions can be approximated in Eq. (4.8) by a uniform value that we call p script  (r). That is, p&#39;(r&#34;) in Eq. (4.8) can be replaced by p script  (r)θ i  (r&#34;), and we find 
     
         P&#39;(r)=fp.sub.script (r),                                   (4.9) 
    
     
         where 
    
     
         f=∫Δ(r-r&#34;)θ.sup.1 (r&#34;)dr&#34;                 (4.10) 
    
     is the macroscopic fill fraction of inclusions, here for simplicity assumed to be essentially uniform throughout the medium. To find an expression for P script  (r), we return to Eq. (4.6). 
     Look first at the term involving T c  : For a point r&#34; in a given inclusion, the integral over points r&#39; in the same inclusion will give zero by virtue of Eq. (4.4b), the assumption of a mesoscopically uniform p script  (r), and the fact that c(r&#34;-r&#39;) is essentially uniform over consider the integral over points r&#39; in different inclusions. Since the range R of c(r&#34;-r&#39;) satisfies R&gt;&gt;b, there are many neighboring inclusions involved. If these inclusions are randomly distributed with respect to the first inclusion, we expect on average no contribution from the integral, by virtue of the assumption of a mesoscopically uniform p script  (r) and Eq. (4.4a). Neglecting the corrections that could result from any correlations in the positions of the inclusions, we set the total contribution from the integral involving T c  equal to zero. 
     Next, note that in Eq. (4.6) E c  (r&#34;) is already a macroscopic field, obtained from the spatial average over a range R of the mesoscopic field e(r) +4πp&#39;(r)3ε h . Thus, its variation over distances of order R can be expected to be small, and we can set E c  (r&#34;)≃E c  (r) in Eq. (4.6) for points r&#34; within R of r. It is then clear that the mesoscopically uniform polarization p script  (r) appearing in Eq. (4.9) should be taken as 
     
         4π.sub.script (r)=3ε.sup.h βE.sup.cl (r).  (4.11) 
    
     Using Eq. (4.11) in (4.9), and recalling that here ##EQU15## (Eq. 4.2 with P NL  =0), we find that 
     
         4πP&#39;(r)=3ε.sup.h βf(1-βf).sup.-1 E(r).(4.13) 
    
     Finally, using Eq. (2.6) in Eq. (2.4) we have, neglecting nonlinear effects, 
     
         p(r)=χ.sup.h e(r)+p&#39;(r) 
    
     
         P(r)=χ.sup.h E(r)+P&#39;(r),                               (4.14) 
    
     where the second of Eqs. (4.14) comes from spatial averaging the first. The total displacement D(r) is thus given by 
     
         D(r)=E(r)+4πP(r)=ε.sup.h E(r)+4πP&#39;(r)≡εE(r),(4.15) 
    
     where the last of Eqs. (4.15) defines the effective-medium dielectric constant ε. Using Eq. (4.13) in Eq. (4.15), we find that ε satisfies ##EQU16## which is the usual Maxwell Garnett result. 
     We close this section by deducing some relations that will prove useful in the following section. Using Eqs. (4.13) and (4.16) in Eq. (4.12) we find that ##EQU17## Next, we combine Eqs. (4.5), (4.11), and (4.17) to determine the mesoscopic electric field e(r) inside an inclusion; the result is ##EQU18## if r lies in an inclusion. The mesocopic field e(r) at points outside an inclusion is given below in Eq. (6.4). Finally, the dipole moment μ of an inclusion associated with the polarization p script  (r) of Eq. (4.11) is given by ##EQU19## 
     NONLINEARITY IN THE INCLUSIONS 
     We now turn to the nonlinear problem, where p NL  (r)≠0. Combining Eqs. (2.4) and (2.6), we find that the total mesoscopic polarization can be expressed as 
     
         p(r)-χ.sup.h e(r)+p&#39;(r)+p.sup.NL (r),                  (5.1) 
    
     (cf. Eq. 4.14). Spatial averaging of this result gives 
     
         P(r)=χ.sup.h E(r)+P&#39;(r)+P.sup.NL (r),                  (5.2) 
    
     and introducing the definition of D(r) we have 
     
         D(r)≡E(r)+4πP(r)=ε.sup.h E(r)+4πP&#39;(r)+4πP.sup.NL (r).(5.3) 
    
     For a medium that is weakly nonlinear, we expect that, when P&#39; (r) and P NL  (r) are written to third order in the electric field strength E(r), we will have a relation of the form 
     
         D(r)=εE(r)+4πA[E(r)·E*(r)]E(r)+πB[E(r)·E(r)]E*(r).                                                    (5.4) 
    
     where the tensor nature of the third-order terms in Eq. (5.4) follows from symmetry arguments if, at the macroscopic level, the medium is isotropic. We neglect nonlinear processes such as third-harmonic generation, and concern ourselves only with the nonlinear modification of the propagation of light at its incident frequency. The two parameters A and B are in general independent, although for electronic nonlinearities we must have A-B→0 as ω→0. The determination of A and B in terms of the corresponding nonlinear coefficients of the host and inclusions is the problem to which we now turn. 
     In the present section, we treat the problem in which only the inclusion material responds nonlinearly to the optical field. Assuming for simplicity that the inclusion material is itself isotropic and uniform, we then have a nonlinear polarization of the form 
     
         P.sup.NL (r)=θ.sup.i (r)p.sup.NL,i (r),              (5.5) 
    
     
         where 
    
     
         P.sup.NL,i (r)=A.sup.i [e(r)·e*(r)]e(r)+1/2B.sup.i [e(r)·e(r)]e*(r),                                (5.6) 
    
     and where A i  and B i  are the nonlinear coefficients of the inclusion material. Since we are looking only for the lowest-order macroscopic nonlinearity (Eq. 5.4), it suffices to estimate e(r) in Eq. (5.6) from the results of a calculation that neglects the nonlinearity itself. That is, we use Eq. (4.18) in Eq. (5.6) to obtain ##EQU20## Next, we average Eq. (5.5) (cf. Eqs. 2.12 and 5.7) to find P NL  (r). Using the fact that E(r) can be assumed to vary little over a range R, and using expression (4.10) for f, we find that ##EQU21## 
     Referring back to Eqs. (5.2)-(5.4), we see that, in order to determine D(r) and subsequently to find expressions for A and B, our remaining task is to find P&#39;(r). This result cannot be taken simply from the linear calculation of section 4, because p NL  (r) makes a contribution to e(r) and thus to p&#39;(r) to lowest nonvanishing order in the nonlinearity. From Eq. (4.5) we have 
     
         4πp&#39;(r&#34;)=(ε.sup.i -ε.sup.h)θ.sup.i (r&#34;)e(r&#34;),(5.9) 
    
     but we must now use the full set of Eqs. (4.1) and (4.2), with p s  (r&#34;)=p&#39;(r)+p NL  (r), P s  (r)=P&#39;(r)+P NL  (r). Equation (5.9) thereby becomes ##EQU22## where β is given by Eq. (4.7), and where the cavity field now contains a nonlinear contribution and is given by ##EQU23## It follows from the discussion after Eq. (4. 8) that p&#39;(r&#34;) can be approximated within the inclusion by a mesoscopically uniform value, to be denoted by p script  (r), at points r&#34; within R of r; then as in section 4, 
     
         P&#39;(r)=fp.sub.script (r).                                   (5.12) 
    
     Note that, from Eq. (5.7) and the fact that E(r) varies little over a range of R p NL ,i (r&#34;) has also been taken to be uniform in this sense, p NL ,i (r&#34;)=θ i  (r&#34;)p NL   script  (r). As discussed after Eq. (4.8.) we see that again we may neglect contributions from the integral involving T c , and we find (instead of Eqs. 4.11, 4.12) that ##EQU24## E c  (r) is given by Eq. (5.11). Combining Eqs. (5.7-5.8, 5.11-5.13), we can solve for P&#39;(r); in fact, it is clear from Eq. (5.3) that we need the sum of P&#39;(r) and P NL  -(r), for which we find ##EQU25## where P NL  (r) is given by Eq. (5.8) and c is given by Eq. (4.16). Using Eqs. (5.8) and (5.14) in Eq. (5.3) and comparing with Eq. (5.4), we identify ##EQU26## as the nonlinear coefficients of the effective medium. Note that the &#34;local field correction factor,&#34; in this case (ε+2ε h )/(ε i  +2ε h ), appears in fourth-order in Eqs. (5. 15) ; three powers enter because the nonlinearity is cubic in the field (Eqs. 5.6 and 5.8), while the fourth appears because the material also responds linearly to the field generated by the nonlinear polarization (Eqs. 5.13 and 5.14). 
     The results given by Eqs. (5.15) are illustrated graphically in FIG. 3. There the vertical axis can be taken to represent either the nonlinear coefficient A of the composite normalized to the nonlinear coefficient A i  of the inclusion material or the value of B for the composite normalized to the nonlinear coefficient B i  of the inclusion material. The horizontal axis gives the fill fraction f of the nonlinear material. For illustrative purposes we have plotted these curves over the entire range 0≦f≦1, even though inter-inclusion correlations not included in this theory can in general be expected to modify significantly the results if we do not have f&lt;&lt;1. For the case in which the linear dielectric constants of the host and inclusion materials are equal, we see that the nonlinear coefficients of the composite are simply equal to those of the inclusion multiplied by the fill fraction f. Thus there is no enhancement in the nonlinear coefficients and of the nonlinear susceptibility. More generally, we see that the nonlinearity of the composite material increases with f at a rate that is either more or less rapid than linear, depending upon the ratio of the linear dielectric constants. This result can be appreciated in that for an inhomogeneous material having the illustrated topology, the electric field will tend to become concentrated in regions of lower dielectric constant. For ε i  &lt;ε h , the electric field within the nonlinear component (the inclusion material) will be larger than the spatially averaged electric field, thus enhancing the effective nonlinearity of the composite material. 
     Note also that Eqs. (5.15) predict that, for the case in which only the inclusion material is nonlinear, the ratio of nonlinear coefficient A/B for the composite is equal to the ratio A i  /B i  of the inclusion material. In the next section, we shall see that in the opposite limiting case in which only the host is nonlinear, the ratio A/B for the composite is not necessarily equal to the ratio A h  /B h  for the host. 
     NONLINEARITY IN THE HOST 
     We now consider the more complicated situation in which the host material is nonlinear with response coefficients A h  and B h , but the inclusion material is linear. Equations (5.1)-(5.4) are then still valid, but instead of Eqs. (5.5) and (5.6) we have 
     
         p.sup.NL (r)=θ.sup.h (r) p.sup.NL,h (r),             (6.1) 
    
     where θ h  (r)=1 if the point r is in the host material and 0 if r is in an inclusion, and where 
     
         p.sup.NL,h (r)=A.sup.h [e(r)·e*(r)]e(r)+1/2B.sup.h [e(r)·e(r)]e*(r).                                (6.2) 
    
     As in section 5, it suffices to estimate e(r) in Eq. (6.2) from the results of a calculation (section 4) performed neglecting the nonlinearity itself. In the linear limit p s  (r)=p&#39;(r), where p&#39; (r) vanishes within the host material, so for use in Eq. (6.2) the field e(r) of Eq. (4.1) reduces to 
     
         e(r)=E.sup.c (r)+∫T.sup.c (r-r&#39;)·p&#39;(r&#39;)dr&#39;.  (6.3) 
    
     Now consider a particular inclusion, centered for simplicity at the origin. In its neighborhood a nonlinear polarization will be induced according to Eqs. (6.1)-(6.3). In evaluating the contribution from the T c  term in Eq. (6.3), we may neglect the contributions from all the other inclusions except the one at the origin, following the discussion in the paragraph following Eq. (4.10). So outside our particular inclusion we have 
     
         e(r)=[0+a.sup.3 βε.sup.h T(r)]·E.sup.c,(6.4) 
    
     where we have used Eq. (4.19) for the (linear) dipole moment of the inclusion, and Eq. (4.4c) for the integral over the inclusion. We have assumed we are at distances r&lt;b&lt;&lt;R from our particular inclusion, so T c  has been approximated as T 0  ; we shall see shortly that it is only at such distances that we get a significant contribution from the T term in Eq. (6.4) to Eq. (6.2). The fact that r&lt;&lt;R has allowed us to treat E c  (r) as essentially uniform E c  ≡E c  (r0). 
     Writing the Cartesian components of Eq. (6.4) as 
     
         e.sub.i (r)=[δ.sub.ij +r.sup.-3 a.sup.3 βt.sub.ij ]E.sup.c.sub.j,(6.5) 
    
     where repeated indices are to be summed over, and where 
     
         t.sub.ij ≡(3n.sub.i n.sub.j -δ.sub.ij),        (6.6) 
    
     with n≡r/r, we can use Eqs. (6.2) and (6.5) in Eq. (6.1) to determine p NL  (r) in the neighborhood of our particular inclusion. We find 
     
         p.sup.NL (r)=p.sup.NLu (r)+p.sup.NLd (r),                  (6.7) 
    
     
         where 
    
     
         p.sup.NLu (r)=A.sup.h [E.sup.c ·E.sup.c* ]E.sup.c +1/2B.sup.h E.sup.c ·E.sup.c ]E.sup.c*,                      (6.8) 
    
     
         p.sup.NLd.sub.m (r)=A.sup.h D.sup.βββ*.sub.mljk (r)E.sup.c.sub.l E.sup.c.sub.j E.sup.c*.sub.k +1/2B.sup.h D.sup.β*ββ.sub.mljk (r)E.sup.c*.sub.l E.sup.c.sub.j E.sup.c.sub.k,                                            (6.9) 
    
     where D.sup.βββ* mljk  and D.sup.β*ββ mljk  are specified by Eq. (II.11) of Appendix II, and the discussion following Eq. (II.12). Note that p NLu  (r) is uniform in the neighborhood of our particular inclusion, while p NLd  (r) is a nonlinear polarization that &#34;dresses&#34; the inclusion, and is nonzero only close to it. Of course, there are term p NLd  (r) that cancel p NLu  (r) when r is inside the inclusion (see Eq. II. 11 of Appendix II), since p NL  (r) is zero there. Now we put Eq. (6.7) into the third of Eq. (2.12) to determine P NL  (e,rar/r/ ). We treat E c  as uniform over the range R of Δ(r) (recall the discussion preceding Eq. 4.11); adding up the contributions from the D terms inclusion-by-inclusion, and neglecting the variation of over the range of the D&#39;s, we find that ##EQU27## where N script  ≡f/(4πa/3) is the number of inclusions per unit volume. The integrals in Eq. (6.10) receive contributions only from terms in the D&#39;s that drop off as r -6  and r -9  for r&lt;a; the integrals are evaluated in Appendix II and are given by Eqs. (II.14). Substituting those expressions into Eq. (6.10), and using Eq. (4.17) for E c , we find that ##EQU28## where β is given by Eq. (4.7) and ε by Eq. (4.16). The terms A h  (1-f) and B h  (1-f) describe the nonlinear polarization that would result if all of the host material (fill fraction (1-f)) responded nonlinearly only to the cavity field; the other terms result from the &#34;dressings of nonlinear polarization&#34; induced above and beyond this in the neighborhood of each inclusion by its own dipole field. 
     To complete our expression (5.3) for the displacement field D(r), and subsequently to identify the effective-medium nonlinear response coefficients A and B (Eq. 5.4), we must now find P&#39; (r). Since p&#39; (r) is nonzero only in the inclusions (see Eq. 2.6), we return to the consideration of one particular inclusion as discussed at the start of this section. For points r at r&lt;a we have ##EQU29## where we have used Eqs. (2.6), (2.8), and (4.1), as well as the fact that p NL  (r)=0 inside the inclusion; E c  of course is given by (see Eq. 4.2) ##EQU30## where we again neglect the variation of E c  over our particular inclusion. The integral in Eq. (6.12) involves p s  (r)=p NL  (r)+p&#39;(r). Following the arguments given after Eq. (4.10), the contribution from p&#39;(r) is negligible and, arguing similarly with p NL  (r) we can neglect the contributions from the dressings of nonlinear polarizations (see the discussion after Eq. 6.11) surrounding all inclusions other than the particular one under consideration That is, the p s  (r&#39;) in Eq. (6.12) can be replaced by the p NL  (r&#39;) given by Eq. (6.7). Doing this, and integrating Eq. (6.12) over our particular inclusion we find ##EQU31## is the dipole moment associated with the p&#39;(r) of our particular inclusion, and where 
     
         E.sub.script ≡∫T(r&#39;)·p.sup.NL (r&#39;)dr&#39;, (6.16) 
    
     with p NL  (r&#39;) given by Eq. (6.7). In arriving at Eqs. (6. 14) and (6. 16) we have omitted the cut-off function c(r) from T c . For p NLu  it is not needed, as long as we integrate Eq. (6.16) over solid angle first, since p NLu  is uniform; and p NLd  is confined to the neighborhood of the inclusion, so a cut-off function of range R&gt;&gt;b is not required. In arriving at the form of Eq. (6.16) we have also been able to use Eq. (4.4c), since p NL  (r&#39;) vanishes for points r&#39; note also that T(-r)=T(r&#39;). As usual, we have neglected the variation of macroscopic fields over distances on the order of a &lt;&lt;R. Putting Eqs. (3.5) and (6.7) in Eq. (6.16) we may evaluate E script  ; we find ##EQU32## where the definitions of the integrals appearing in Eq. (6.17) are given in Appendix II; the values of the integrals are also worked out there (Eq. II.13). With these results, we find ##EQU33## 
     We now determine P&#39;(r) by averaging p&#39;(r). Neglecting the variation of Δ(r) over distances on the order a, we find that the integral (2.12) reduces to essentially a sum of the dipole moments μ&#39;. Using Eq. (6.14) and neglecting as usual the variation of macroscopic fields over distances of the order of R, we find 
     
         4πP&#39;(r)=3ε.sup.h βf[E.sup.c (r)+E.sub.script (r)],(6.19) 
    
     (contrast with Eq. 4.11), were now E script  (r&#39;) is given as a function of r by Eq. (6.18). Putting Eq. (6.13) into Eq. (6.19) we find that we may write ##EQU34## where ε is given by Eq. (4.16), and we have used Eq. (4.17) for β; using Eq. (6.20) in Eq. (5.3) we can then write 
     
         D(r)=εE(r)+D.sup.NL (r),                           (6.21) 
    
     
         where ##EQU35## Substituting Eqs. (6.11) and (6.18) into (6.22) and comparing with Eq. (5.4), we may identify the nonlinear response coefficients of the effective medium. We find ##EQU36## As in the results of the preceding section, there are four &#34;local field correction factors&#34; which appear in the expressions above. They are different here than in the preceding section, because they apply to the local field in the host (see Eq. 4.17) rather than in the inclusion (see Eq. 4.18). Beyond that, the results for the nonlinearity in the host are much more complicated because, as is clear from comparing the derivations in sections 5 and 6, the mesoscopic fields are much more complicated in the host. In particular, while if the nonlinearity is in the inclusions we have A/B=A.sup.i /B.sup.i, we do not have A/B=A.sup.h /B.sup.h if the nonlinearity is in the host. Plots of different cases which apply the results of Eqs. (6.23) are shown in FIGS. 4 and 5, (a), (b), (c). 
    
     In FIG. 4 we plot the enhancement in the quantity A+(1/2) B [i.e. we plot (A+(1/2)B)/(A h  +(1/2)B h )] as a function of the fill fraction f of inclusion material for several values of the ratio ε i  /ε h  of linear dielectric constants. We have the quantity A+(1/2)B because it is proportional to the nonlinear refractive index experienced by linear polarized light; as mentioned above, in general the enhancement in A is different from the enhancement in B. Moreover, as can be verified by careful examination of Eqs. (6.23). The enhancement in A+(1/2)B for a given value of f does not depend upon the ratio B h  /A h  of nonlinear coefficients of the host material; conversely, the enhancement in A or B separately does depend on the ratio B h  /A h . For the case in which the inclusion and host have the same value of the dielectric constant, we see that A+(1/2)B decreases linearly to zero as f increases from zero to one. The origin of this behavior is simply that the amount of nonlinear material decreases as the fraction of inclusion material increases. Nonetheless, for ε i  /ε h  greater than one and preferably greater than approximately two, we find that the quantity A+(1/2)B increases as the fill fraction of (linear) inclusion material increases. This effect occurs because the presence of the inclusion material modifies the electric field distribution within the composite material in such a manner that the spatially averaged cube of the electric field within the host material is significantly increased. 
     Some of the tensor properties of the nonlinear response of the composite material are shown in FIGS. 5a, b and c. In each graph, the ratio B/A for the composite is plotted as a function of the fill fraction f for several values of the ratio ε i  /ε h  of linear dielectric constants. FIG. 5a, FIG. 5b and FIG. 5c refer respectively to the cases B h  =0, B h  /A h  =1, and B h  /A h  =6, and correspond physically to a nonlinear response dominated by electrostriction (B h  =0), electronic response in the low-frequency limit (B h  /A h  =1), and molecular orientation (B h  /A h  =6). From part (a) we see that, even when B h  vanishes, the composite can possess nonlinear coefficients A and B that are comparable in size. From part (b) we see that, for the special case B h  /A h  =1, the ratio of nonlinear coefficients of the composite is equal to that of the host for any value of f. 
     SUMMARY 
     In summary, we have generalized the Maxwell Garnett Theory of the optical response of composite materials by allowing either or both constituents of the material (the host component and the inclusion component) to possess a third-order nonlinear susceptibility. Equations (5.15) and (6.23) respectively give the key results of our calculation for cases in which only the inclusion material and in which only the host material is nonlinear. When both components respond nonlinearly, the effective values of the nonlinear susceptibilities A and B are obtained by summing the two contributions given by Eqs. (5.15) and (6.23). In these equations, A and B are defined by Eq. (5.4), A i  and B i  by Eq. (5.6), A h  and B h  by Eq. (6.2), ε by Eq. (4.16) and β by Eq. (4.7). Our treatment takes full account of the tensor nature of the nonlinear interaction under the assumptions that each component is optically isotropic and that the composite is macroscopically isotropic. For the case in which only the inclusion material is nonlinear, our results are consistent with those set forth in the articles by D. Rich et al., G. S. Agarwal et al., and J. W. Haus et al., referenced above. For the case in which the host material is nonlinear, the nonlinear susceptibilities A and B for the composite can be considerably larger than those of the host material itself; moreover, the ratio B/A for the composite can be very different from that of the host. We note that in all of the examples presented we have assumed that both constituents are lossless and thus that the parameters ε i , ε h , A i , B i , A h , and B h  are all real. Nonetheless, the formulas presented here are correct even in the more general case where those parameters are complex. 
     SYSTEMS DEVICES AND METHODS 
     Systems in which nonlinear interactions occur in accordance with the invention utilize the material system discussed above and means for applying optical radiation to establish a field in the material system (the body of composite material). The field is applied by propagating a laser beam or beams so that they are incident on the body, and preferably are transmitted therethrough. Configurations of the system for different nonlinear optical interactions are described in Section 1.2 of the text by Robert W. Boyd referenced above. 
     For second harmonic generation a beam of frequency ω propagates through the body which is a composite material system having second-order nonlinear susceptibility, χ.sup.(2). For some frequency generation two beams of different frequencies ω 1  and ω 2  are applied to the body having the χ.sup.(2) nonlinear susceptibility. Similarly, two beams of different frequency produce an output beam of difference frequency when the intensity of the ω 1  beam exceeds that of the χ 2  beam. The sum frequencies are provided when the ω 1  and ω 2  beams are of approximately equal intensity. Optical parametric oscillation is obtained by containing the body of nonlinear material in a cavity and pumping the body with a laser beam. The refractive index changes as a function of the intensity of the incident beam and provides self-focusing (a lens) in the body of composite material. 
     Third order material, χ.sup.(3), is used for third harmonic generation and also when a plurality of beams are applied so as to obtain sum or difference frequency generation depending upon the relative intensities of the incident beams. Photonic devices, for example an optically bistable device, is obtained by utilizing the body as a saturable absorber. Other photonic and electro-optic devices and systems are obtained by utilizing the body within a control electric field as discussed in Chapter 10 of the above-referenced text by Robert W. Boyd. In all cases the efficiency and the nonlinear interactions in the body are enhanced because of the enhanced susceptibility obtained by virtue of the host being of material having nonlinear optical susceptibility (a nonlinear response to the amplitude (intensity) of the incident laser beam or beams). 
     MATERIAL SYSTEMS 
     Material systems embodying the invention have the topology of inclusion particles imbedded in the host. Alternatively, the systems may have two interdispersed components. Where the inclusion particles are provided by a matrix containing the nonlinear component, the nonlinear component then acts as the host in the case where the inclusion particles are imbedded in the host (FIG. 1). 
     In the following discussion, systems which provide enhanced third-order susceptibility (χ.sup.(3)) will first be discussed, then systems which provide enhanced second-order nonlinear susceptibility (χ.sup.(2)) will be discussed. The first system consists of particles suspended in carbon disulfide (CS 2 ). The particles may be titania (TiO 2 ). Titania is presently preferred because of its high refractive index (2.8) and because it may be precipitated as particles in the preferred particle size range (the order of 50 Angstroms). Carbon disulfide is preferable because of its low refractive index (1.6) and because it is a well accepted nonlinear optical material. The host and particles thus differ significantly in their respective refractive indices. Titania has nonlinear susceptibility characteristics. 
     The particles may alternatively be diamond particles. Semiconductor microcrystallites may also be used. It is desirable to use semiconductors, such as cadmium sulfide at photon energies (wavelengths of the incident beam) well below the band gap energy of the semiconductor. For example, cadmium sulfide microcrystallites are desirably used when the beam is in the infrared frequency range because band gap energy corresponds to a wavelength of approximately 700 Angstrom units and the near infrared is from approximately 1.3 to 1.55 micrometers in wavelength. There the semiconductor is transparent and lossless to the incident irradiation. 
     Material systems which emphasize the third-order nonlinearities or susceptibilities utilize hosts which are isotropic. Material systems which emphasize the second-order nonlinear response or susceptibility are desirably anisotropic. A preferred host material for systems operating in accordance with third-order nonlinear optical interactions may be plastics and polymers. Conjugated polymers having small linear refractive indices are preferred. An especially preferred polymer is polydiacetylene. The chemistry of such polymers and various types of conjugated polymers are described in the text by P. N. Prasad and D. J. Williams, Introduction to Nonlinear Optical Effects in Molecules and Polymers, John Wiley &amp; Sons, Inc., New York (1991). A listing of such polymers is found in Table 10.02 in that text. In addition to polydiacetylenes, other polymers which may be used as host components of the material system are poly-ρ-phenylene vinylenes, polyacetylene, polythiophene, polyimides, phthalocyaine, polysilanes and ferocene. 
     In addition to carbon disulfide, other liquids (at room temperature) which exhibit third-order nonlinear response are listed in Table 10.1 of the above-referenced text. They include benzene, allo-ocimine, trans-retiliol, nitrobenzene, iodobenzene, pyridine, and thiophene. 
     The particles may be suspended in the host, especially in the case of plastics while the plastic is in its liquid state and the particles are held in place in the plastic when it cures to solid state. The composite may be in the form of a colloid. Then the particles are suspended in the nonlinear host medium. It is possible to use metal particles (gold, silver or platinum) in a nonlinear host medium in order to provide nonlinearity. Then absorption in the particle of the illumination must be tolerated. 
     The interdispersed topology may be implemented using porous glass or glass which is formed from sol gel media. One porous glass is sold by Corning Glass Company of Corning, N.Y. under the trademark Vycor. The host is introduced, as by capillary action, into the pores and is interdispersed with the material of the pores which then provide the effect of the particles. 
     Second-order nonlinear systems preferably consist of poled polymers which are cured in the presence of an electric field. Such polymers are discussed in the above-referenced text by Prasad and Williams and their chemistry and physics is described in Section 4.2 of the text. The polymers may be constructed as multi-layers with the particles distributed between the layers. Formation of multi-layers from plastic and other polymeric materials is known in the art and has been used to provide Langmuir-Blodgett films. See Section 4.4 of the Prasad and Williams text. Poled polymers are described in Section 7.4 of the text. Suitable material for forming poled polymers are p-(dimethylamino)-p*-nitrostilbene, polystyrene, PMMA, parachlorosubstituted polymer, copolymers with methyl methacrylate or side chain substituted methacrylates, acrylic polymers, diaminodicyanoquinodimethanes polycarbonate with thiapyrylium disalt, and polyethyleneoxide and p-nitroaniline dissolved in poly(ε-caprolactone). The structural formulas and methods for producing these materials are found in the Prasad and Williams text or references cited therein. 
     From the foregoing description, it will be apparent that there has been provided significant improvements in nonlinear optics and particularly in systems, devices and methods of processing optical radiation (beams and signals) which exhibit enhanced nonlinear response to the radiation. Variations and modifications in the herein described systems, devices and methods, within the scope of the invention, will undoubtedly suggest themselves to those skilled in the art. Accordingly, the foregoing description should be taken as illustrative and not in a limiting sense. 
     APPENDIX I 
     In this appendix we derive a relation between the electric field e(r) satisfying the Maxwell equations 
     
         ∇·[εe(r)=-4π∇·p(r) 
    
     
         ∇·b(r)=0 
    
     
         ∇×e(r)-iωb(r)=0 
    
     
         ∇×b(r)+iωεe(r)=-4πiωp(r),(I.1) 
    
     and its macroscopic average E(r), 
     
         E(r)=∫Δ(r-r&#39;)e(r&#39;)dr&#39;.                          (I.2) 
    
     Equations (I.1) are identical to Eqs. (2.7), with ε h  and p s  (r) replaced by ε and p(r) respectively to simplify the notation; ε is taken to be uniform and, as before, ω≡ω/c. For most of this appendix (up to and including Eq. I.43) we need assume only that Δ(r) is a spherically symmetric function, Δ(r)=Δ(r), where r=|r|, that vanishes sufficiently rapidly as r→∞. 
     A relation similar to the one we find here was considered earlier [14], where we took ε=1 and assumed that p(r) was a sum of Dirac delta functions. Here we are interested in ε≠1, and in a p(r) that is continuous except for stepwise discontinuities at the surfaces of our inclusions; in the usual way p(r) is assumed to vanish as r→∞, and we are interested in solutions of Eqs. (I.1) for which the particular component satisfies the usual outgoing radiation condition [17]. 
     We begin with two forms of the solution to Eqs. (I.1). The first is obtained by noting that those equations can be written as ##EQU37## where k≡√εω is taken to have Im k≧0, Re k&gt;0 if Im k=0. Thus the solutions for εe(r), √εb(r) may be recovered, upon replacement of ω by k, from the solutions in the ε=1 limit. The latter solutions are well-known [18], and so we identify from the solution to Eq. (I.1) for e(r) as ##EQU38## where U is the unit dyadic; in component form ##EQU39## where the subscripts label Cartesian components. The notation in Eq. (I.4) indicates that a small sphere of radius η, centered at r&#39;=r, is to be excluded from the r&#39; integration, and then at the end of the calculation the radius η of the excluded sphere is to be allowed to approach zero. This &#34;excluded sphere&#34; is necessary to make Eq. (I.4) even unambiguous [19], since F(r) diverges (as r -3 ) as r→0. A convenient shorthand is to rewrite Eq. (I.4) as ##EQU40## The expression (I.4) or (I.7) is valid at all points r except those (a set of measure zero) where p(r) is changing discontinuously; at such points e(r) is itself varying discontinuously in a steplike manner, as described by Eq. (I.4) or (I.7) as r moves from one side of the surface to the other. Finally, we note that we have omitted any homogeneous solution of Eqs. (I.1) in writing down Eq. (I.4). We continue to omit such homogeneous solutions until later in this appendix (see Eq. I.44). 
     A second form of the particular solution to Eq. (I.1) may be obtained by taking the curl of the third of those equations, using the identity curl curl=grad div-∇ 2 , , and the first of Eqs. (I.1) to find 
     
         [∇.sup.2 +k.sup.2 ]e(r)=-4πq(r),               (I.9) 
    
     
         where 
    
     
         q(r=ω.sup.2 p(r)+ε.sup.-1∇(∇· p(r).(I.10) 
    
     The particular solution to Eq. (I.9) is well known [20], ##EQU41## Strictly speaking, an excluded sphere centered at r&#39;=r should be specified in Eq. (I.11), since the Green function |/r-r&#39;| -1  exp(ik|r-r&#39;|) is undefined as r&#39;→r [16]. But because the divergence is in practice masked by the volume element dr, the omission of the excluded volume here usually does not lead to difficulties. We proceed now by using Eq. (I.11) to find an expression for E(r). 
     Using Eq. (I.11) in Eq. (I.2), we find 
     
         E(r)=∫I(k;r-r&#39;)q(r&#39;)dr&#39;.                              (I.12) 
    
     
         where ##EQU42## Using the Green function expansion [21] ##EQU43## where r.sub.&lt;  and R.sub.&gt;  are respectively the lesser and greater of r and r&#39;, the Y.sub.lm are the spherical harmonics, and j.sub.l and h.sup.(1).sub.l are respectively the l.sup.th order spherical Bessel and Hankel functions of the first kind, we find ##EQU44## which shows that I(k;r) depends only on k and r=|r|. Writing the first integral as the difference between an integral over all r&#39; and one between r and ∞, we use h.sub.o.sup.(1) =(ikr).sup.-1 exp(ikr) to find ##EQU45## with n a unit vector in an arbitrary direction, and ##EQU46## The second form of Eq. (I-18) is obtained by using h.sub.o.sup.(1) (x)=j.sub.o (x)+in.sub.o (x), where n.sub.o (x) is the zeroth-order spherical Bessel function of the second kind; it demonstrates that I(k;r) is purely real. 
    
     Inserting Eq. (I-16) into Eq. (I.12) we find that ##EQU47## It is now convenient, for a given field point r, to write 
     
         p(r&#39;)=p.sup.&gt; (r&#39;)+p.sup.&lt; (r&#39;),                           (I.21) 
    
     where ##EQU48## and where for the moment η is an arbitrary positive number. We then put ##EQU49## and likewise for q &lt;  (r&#39;), e &lt;  (r), E &lt;  (r); here ∇&#39; ≡∂/∂r&#39;. Note that q &lt;  (r&#39;) and q &gt;  (r&#39;) are both singular |r-r&#39;|=η, although their sum is not. Nonetheless, since p &gt;  (r&#39;) vanishes at r&#39;=r, we may insert the first of Eq. (I.23) into the first integral in the third of Eqs. (I.23) and perform two partial integrations. The result is 
     
         E.sup.&gt; (r)=Δ(k)∫F(r-r&#39;)·p.sup.&gt; (r&#39;)dr&#39;+E.sup.(1) (r)+E.sup.(2) (r)                                         (I.24) 
    
     
         where 
    
     
         E.sup.(1) (r)=ε.sup.-1 ∫I(k;r-r)∇&#39;(∇&#39;·p.sup.&gt; (r&#39;))dr&#39; 
    
     
         E.sup.(2) (r)=ω.sup.2 ∫I(k;r-r&#39;)p.sup.&gt; (r&#39;)dr&#39; (I.25) 
    
     and the F in Eq. (I. 24) is given by Eq. (I. 5). We next partially integrate the express E.sup.(1) (r) twice, again using the fact that p&gt;(r&#39;), and ∇&#39;·p &gt;  (r&#39;), vanish at r&#39;=r. The first partial integration yields 
     
         E.sup.(1) (r)=ε.sup.-1 ∫∇I(k;r-r&#39;)(∇&#39;·p.sup.&gt; (r&#39;))dr&#39;(I.26) 
    
     using γ&#39;I(k;r-r&#39;)=-∇I(k;r-r&#39;). From Eq. (I.18) we find ##EQU50## with n&#39; o  (x)≡dn o  (x)/dx and j&#39; o  ≡dj o  (x)/dx, and where we have used ##EQU51## (see Eqs. 1.19). The second partial integration gives 
     
         E.sup.(1) (r)=ε.sup.-1 ∫(p.sup.&gt; (r&#39;)·∇)∇(k;r-r&#39;).              (I.30) 
    
     To evaluate this expression we need the components ##EQU52## with n&#34; o  (x)=dn&#39; o  (x)/dx, etc.; we have used Eq. (I.28) and the Wronskian relation j o  (x)n&#39; o  (x)-j&#39; o  (x)n o  (x0=x -2 , Using Eq. (I-32) in Eq. (I-31) we write ##EQU53## Using Eq. (I.35) in Eq. (I-24) along with the second of Eqs. (I.25), we find ##EQU54## Finally using Eqs. (I.18) and (I.35) in Eq. (I.39) we can write ##EQU55## 
     We now return to Eq. (I.22) and consider the effect of taking the limit η→0. For certain fields this must be done with some care. For example, it is not true that e &lt;  (r)→0 as η→0, because the field inside a polarized sphere centered at r will not vanish even as the radius of that sphere becomes vanishingly small. However, it is true that E &lt;  (r)→0 as η0, since the fraction due to the sphere of the volume over which e &lt;  (r) is integrated to yield E &lt;  (r) vanishes as η→0, and of course the field outside the small sphere of polarization, but due to it, vanishes as η→0; thus as η→0 we have E &gt;  (r)→E(r). We also have P &gt;  (r)→P(r) as η→0, and ##EQU56## which follows from Eq. (I.4). So in the limit η→0 Eq. (I.38) yields ##EQU57## where, following the convention of Eqs. (I.7) and (I.8), we have put ##EQU58## Now recall that in writing Eqs. (I.4), (I.7), and (I.11) we have neglected a homogeneous solution of Eqs. (I.1). Such a homogeneous solution e o  (r) satisfies 
     
         [∇.sup.2 +k.sup.2 ]e.sup.o (r)=0,                 (I.45) 
    
     (cf. Eq. I.9), and will lead to a contribution E o  (r), 
     
         E.sup.o (r)=∫Δ(r-r&#39;)e.sup.o (r&#39;)dr&#39;,            (I.46) 
    
     to the macroscopic field. From Eq. (I.45) we see that e o  (r) will be a sum (or integral) of fields of the form exp(ikn·r), where n is an arbitrary unit vector; so from Eqs. (I.17) and (I.46) we see that 
     
         E.sup.o (r)=Δ(k)e.sup.o (r).                         (I.47) 
    
     Thus, even if the total E(r) includes E o  (r) , and the total e(r) includes e o  (r). Eq. (I.43) will still hold. 
     The result (I.43) of our manipulation is an exact result. It follows directly from the Maxwell Eqs. (I.1), the condition of outgoing radiation implicit in the solutions (I.4) and (I.11), and the fact that the averaging function Δ(r) has spherical symmetry, Δ(r)=Δ(r). We now approximate Eq. (I.43) using the assumption that the range R of Δ(r) satisfies kR &lt;&lt;1. Since Δ 1  (k;r) and Δ 2  (k;r) will also have a range on the order of R, N(k;r) will as well (see Eqs. I.40 and I.41). We can recover an approximation for N(k;r) in the limit kR&lt;&lt;1 by expanding j o  (kr)=(kr) -1  sin(kr) and n o  (kr) =-(kr) -1  cos(kr) in powers of kr to obtain asymptotic series approximations for Δ 1  (k;r) and then using those in Eq. (I.41). For Δ 1  (k;r) we find ##EQU59## and, putting Δ 3  (k;r)≡kr Δ 2  (k;r), we find 
     
         Δ.sub.3 (k;r)-Δ.sub.3 (r)-1/2(kR).sup.2 Δ.sup.II.sub.3 (r)+                                                      (I.50) 
    
     where ##EQU60## Using Eqs. (I.48) and (I.50) in Eqs. (I.41), and expanding the spherical Bessel functions and their derivatives for small kr, we find 
     
         N.sup.(1) (k;r)˜N.sup.(1)(r)+(kR).sup.2 N.sup.(1) II(r)+(I.52) 
    
     where ##EQU61## We are interested only in the lowest order in kR, so we take N(k;r)→N.sup.(1) ·(r). Consistent with this we take (now assuming that Δ(r) is normalized to unity, (Eq. 2.10)), Δ(k)→1. Our exact result (I.43) then simplifies to the approximate expression ##EQU62## is the &#34;cavity field&#34; that is introduced immediately in more heuristic derivations (see section 2), and 
     
         T.sup.c (r)=T.sup.o (r)c(r),                               (I.58) 
    
     where ##EQU63## is the static dipole-dipole coupling tensor in a background medium of dielectric constant ε, cut off for r&lt;η at the origin, and ##EQU64## is a function that cuts off at large r, c(0)=1, c(r)→0, as r→∞, with a range on the order of R. 
     APPENDIX II 
     In this appendix we state some geometrical formulas that are useful in calculating the nonlinear response of the host material in the neighborhood of an inclusion (see section 6). Defining 
     
         r=xx+yy+zz 
    
     
         n≡r/r=x sinθ cosφ+y sinθ sinφ+z cosθ,(II.1) 
    
     where r=|r| and (θ, φ) are the usual angular spherical coordinates, we denote the Cartesian components of n by n 1  : 
     
         n.sub.1 =sinθ cosφ 
    
     
         n.sub.2 =sinθ sinφ 
    
     
         n.sub.3 =cosθ.                                       (II.2) 
    
     For any quantity q=q(θ, φ), we define ##EQU65## where dΩ=sinθdθdφ is an element of solid angle. It is then easy to verify that 
     
         n.sub.i =0 
    
     
         n.sub.i n.sub.j =1/3δ.sub.ij 
    
     
         n.sub.i n.sub.j n.sub.k =0 ##EQU66## where δ.sub.ij is the Kronecker delta (δ.sub.ij =1 if i=j, δ.sub.ij =0 otherwise). All of the results of this appendix are based on Eqs. (I1.4). For example, putting (cf. Eq. 6.6) 
    
     
         t.sub.ij≡ 3n.sub.i n.sub.j -δ.sub.ij 
    
     
         s.sub.ij ≡t.sub.ik t.sub.kj =3n.sub.i n.sub.j +δ.sub.ij,(II.5) 
    
     where throughout this appendix repeated indices are summed over, we easily find ##EQU67## which are used in deriving the results below. 
     We first define some auxiliary quantities that are useful in intermediate steps of the calculations: 
     
         C.sup.αβγ.sub.mljk ≡(δ.sub.ml +σ.sup.α t.sub.ml)(δ.sub.ij +σ.sup.β t.sub.ij)(δ.sub.ik +σ.sup.γ t.sub.ik), 
    
     
         σ.sup.α≡r.sup.-3 a.sup.3 α,        (II.7) 
    
     where the Greek indices label a set of parameters α that are independent of r, θ, and φ; a is also a constant. Further, we define 
     
         A.sup.αβγ.sub.pljk =r.sup.-3 t.sub.pm C.sup.αβγ.sub.mljk 
    
     
         F.sup.αβγ.sub.mljk =C.sup.αβγ.sub.mljk -δ.sub.ml δ.sub.jk.                           (II.8) 
    
     We can then easily determine that ##EQU68## Now our quantities of interest are ##EQU69## where θ(x)=0, 1 for x&lt;0, x&gt;0 respectively. In particular, we seek ##EQU70## where as usual dr=r 2  drdΩ. Using the results (II.9 and II.10), the integrals of B and D in Eq. (II.12) can be easily found. We require two particular cases for the choices of α, β, γ. In the first we want α=β, γ=β*; these terms we denote by B.sup.βββ* pljk  and D.sup.βββ* mljk . In the second we want γ=β, α=β*; these terms we denote by B.sup.β*ββ pljk  and D.sup.β*ββ mljk . Collecting contributions, we find ##EQU71##