A

Measuring two-photon orbital angular momentum entanglement

Abstract

We put forward an approach to estimate the amount of bipartite spatial entanglement of down-converted photon states correlated in orbital angular momentum and the magnitude of the transverse (radial) wave vectors. Both degrees of freedom are properly considered in our framework which only requires azimuthal local linear optical transformations and mode selection analysis with two fiber detectors. The coincidence distributions predicted by our approach give an excellent fit to the distributions measured in a recent experiment aimed to show the very high-dimensional transverse entanglement of twin-photons from a down-conversion source. Our estimate for the Schmidt number is substantially lower but still confirms the presence of high-dimensional entanglement.

pacs:
03.67.Mn, 42.50.Dv, 42.65.Lm

I Introduction

The phenomenon of quantum entanglement whereby distant systems can manifest perfectly random albeit perfectly correlated behavior is now recognized as the essential ingredient to perform tasks which cannot be accomplished with classically correlated systems (1). The presence of entanglement has been traditionally revealed in the violation of Bell-type inequalities (2). However, detecting such a violation does not provide in general a measure of the amount of entanglement. This is particularly significant in systems correlated in multidimensional degrees of freedom (3). Several techniques have been proposed to assess the presence of entanglement for different quantum scenarios. These include state tomography (4); (5); (6), that yields a complete reconstruction of a quantum state but requires many setting measurements, entanglement witnesses (7), which detect some entangled states with considerably less measurements, and the experimental determination of concurrence (8); (9). In Ref. (9), by using two copies of a down-converted two-photon state entangled in polarization and transverse momentum, the measurement of concurrence was achieved with a single, local measurement on one of the photons.

Although bipartite entanglement is well understood, finding experimentally feasible procedures to quantify it for systems correlated in multidimensional degrees of freedom turns out to be quite challenging and relevant. Indeed, access to higher dimensional Hilbert spaces in which information can be encoded and manipulated has recently attracted great interest, with proof-of-principle demonstrations using quantum communication protocols in three-level systems (qutrits), such as entanglement concentration (10), quantum bit commitment (11) and quantum coin-tossing (12). Likewise, a complete characterization of states hyperentangled in polarization, orbital angular momentum and frequency has been experimentally implemented (13).

The aim of this paper is to address the problem of how, by performing a certain set of local linear optical operations affecting one of two multidimensional spatial degrees of freedom in which two-photon states can be entangled, it is possible to obtain an explicit measure of the amount of bipartite transverse entanglement. Specifically, according to the interplay between both spatial degrees of freedom -orbital angular momentum (OAM) and the magnitude of the radial wave vectors- fundamentally different predictions are expected in the subsequent joint photodetection process, via mode-selection analysis with two fiber detectors preceded by azimuthal transformations (acting only on the OAM). The application of our framework is compared with the results of a recent experiment (14).

The paper is organized as follows: Section II presents the general setting of the problem and includes the Schmidt decomposition technique for describing bipartite spatial entanglement. In Section III we reveal the generic features that arise in the photodetection coincidences according to the transverse structure of the two-photon wave function. Section IV illustrates the results found in Section III with an example of a realistic two-photon wave function that yields a full analytical solution to the problem. In Section V a closed-form expression for the Schmidt number is obtained in terms of easily accessible experimental parameters. We then proceed to exploit this Schmidt number in two interesting examples of optical transformations with azimuthal phase plates. These enable us to estimate the amount of spatial entanglement of a two-photon source. Conclusions of the paper are drawn in Section VI. Details of our calculations, together with some useful background material, have been included in two Appendices.

Ii Modal Schmidt Decomposition for Two-Photon States

Two-photon pure quantum states are described in a Hilbert space by a continuous bilinear superposition of spatiotemporal multimode states. Sources of such nonclassical states of light are mostly realised in the process of spontaneous parametric down-conversion (15), where an intense quasimonochromatic laser pump illuminates a crystal endowed with a quadratic nonlinearity producing pairs of photons (idler and signal). Conservation of energy and momentum impose that the state be spectrally and spatially correlated. Here we explore entanglement involving spatial degrees of freedom that depend on the transverse structure of these states. For simplicity, we assume that the down-converted photons are linearly polarized, monochromatic and frequency degenerated. The two-photon state can then be written as , where are the transverse components of the idler and signal wave vectors. Under conditions of paraxial and nearly collinear propagation of the pump, idler and signal photons, the two-photon amplitude is given by (16) . The function represents the transverse profile of the pump, whereas originates from the phase-matching condition in the longitudinal direction and depends on the specific orientation and cut of the nonlinear crystal. Since the arguments of and enforce correlations in different manifolds of the idler and signal wave vector space, it is the global structure of the one dictating the entanglement degree of .

In order to extract the amount of entanglement contained in , one may resort to the Schmidt decomposition that provides the spatial information modes of the two-photon pair. Suppose that the transverse spatial frequency field of the pump beam has a Gaussian profile of the form , where is the width (at the beam waist). The chosen pump profile peaks when its argument vanishes. This imposes that the idler and signal transverse wave vectors should be mostly anticorrelated (). Remarkably, it was shown by Law and Eberly (17) that with such a Gaussian profile the normalised two-photon amplitude can be expressed as

(1)

where and are the normalised polar Schmidt mode functions for the idler and signal photons with topological charge and radial index corresponding to eigenvalues (they satisfy and ) of the reduced density matrices for each photon. Knowledge of yields a direct measure of the degree of transverse entanglement given by the Schmidt number (18) , which is the reciprocal of the purity of the idler and signal density matrices, it is invariant under free propagation and yields an average of the number of relevant spatial modes involved in the decomposition. The larger the value of , the higher the transverse entanglement. For instance, product states correspond to (there is only one nonvanishing eigenvalue equal to ), whereas states with are entangled. A distinguishing feature of decomposition (1) is that it represents in terms of a perfectly correlated discrete basis of paraxial eigenstates of the OAM operator along the direction of light propagation (with corresponding eigenvalues (19); (20), rather than in a continuous plane-wave modal expansion. The precise form of the radial idler and signal eigenmodes and depends on the specific phase-matching function . In particular, when is approximated by a constant there are striking consequences: Eq. (1) becomes a (non-normalisable) superposition in which all OAM eigenstates have an equal weight and the radial dependence becomes just a global factor. It is important to emphasize that to attain decomposition (1) the appropriate choice of the widths, and , for the idler and signal radial eigenmodes, has to be made. If the two-photon amplitude is of the form , then , where is the so-called Schmidt width. For widths different from an additional summation over the radial indexes occurs, and the perfect correlation between the idler and signal radial modes is absent (see Appendix A). Moreover, the fact that the idler and signal mode functions are anticorrelated with respect to their topological charge numbers is a consequence of a more general process: the conservation of OAM, which is transferred from the pump photon (carrying zero OAM for a Gaussian mode) to the down-converted photon pair (21).

Figure 1: (color online). Two-photon coincidence detection configuration. Down-converted idler and signal beams from a nonlinear crystal traverse two optical systems () which include azimuthal phase plates and are coupled into single-mode fibers (SMF) gated by a coincidence circuit (CC).

Iii Azimuthal Transformations on the Two-Photon State

The usefulness of the Schmidt decomposition (1) becomes apparent when analyzing the propagation of two-photon states through optical systems. Each of the intervening modes evolves and transforms independently of the others. It is also clear that the correlation properties displayed by the two-photon amplitude (1) are preserved in the position representation. Therefore, suppose that the idler and signal photon beams, described now by the paraxial two-photon wave function in the transverse position representation (20); (22); (23) , are each transmitted through different linear optical systems that include diffractive (or refractive) azimuthal phase plates (see Fig. 1). The role of the plates in each path is to imprint an azimuth-dependent phase factor on the incoming Schmidt modes. Their (separate) action on the two-photon wave function can be represented by the unitary and radially symmetric impulse response functions . These functions locally transform the spiral harmonic mode content of each photon via , with , so that the resulting output two-photon wave function is

(2)

where the initial perfect correlation in OAM is lost, remaining only that corresponding to the radial modes.

Upon traversing their respective linear optical systems the idler and signal photons, having OAM , and radial indexes , are detected in coincidence. By placing photodetectors at the output ports of suitable arrays of deterministic mode sorter interferometers (24); (25), or computer-generated holograms (10), it is possible to distinguish modes bearing different OAM. In practice, the most straightforward procedure involves projecting into single-mode fibers, where the propagated mode has a fundamental Gaussian profile ().

The probability that the idler and signal photons will be projected into modes and is found to be , with an additional incoherent (and weighted) sum over the radial indexes when taking into account multi-mode detection. This photodection probability is expressed in terms of the spatial overlap of the Schmidt and the fiber modes at the planes where the phase plates are located. Notice that their corresponding widths, and , are not necessarily equal, and thus the orthogonality between these modes when their radial indexes differ does not hold in general. To evaluate , we define , where the two different widths of the radial modes have been specified for clarity. The coincidence probability can then be written as

(3)

All the radial dependence is contained in the functions that modulate the angular impulse response functions . We emphasize that although the radial part of the Schmidt modes does not experience any significant transformation when the idler and signal beams traverse their respective linear optical systems, its proper inclusion in the detection process is essential. This is reflected in the structure of Eq. (3) with the presence of -dependent radial functions , and is a consequence of a non-constant phase-matching function . Had been a constant then the input two-photon wave function (1) would have been represented by a common radial function times a (non-normalisable) maximally entangled superposition of spiral harmonic modes. In this limiting case one derives a fundamentally different prediction for the coincidences: , where now all the radial dependence of the detected modes appears only as a global function.

Iv An Exactly Soluble Model for the Two-Photon Amplitude

To illustrate the previous results, suppose that the phase-matching function is also Gaussian (see Appendix C, where a justification of this model is provided), an approximation implying that photons are generated near the phase-matching region of wave vectors where the down-conversion process occurs most effectively. The normalised two-photon amplitude reads

(4)

where plays the role of an effective width of the phase-matching function, and depends on the nonlinear crystal thickness, although no specific relation is assumed here. At variance with other models, where is often approximated by a constant (), the two-photon amplitude (4) captures the relevant features of the transverse wave-vector correlation between the idler and signal photons, and provides an analytically amenable model that yields explicit formulas for Eq. (3). We show in Appendix A that for the two-photon amplitude (4), the Schmidt mode functions belong to the well-known Laguerre-Gaussian basis (17), or, more generally, to a continuum family of spatial modes generated via metaplectic mappings (rotations on the orbital Poincaré sphere) from the Laguerre-Gaussian modes (20); (26); (27). Furthermore, we also find that the Schmidt eigenvalues are , with . This allows us to express the Schmidt number in the following closed form . Large values of occur when . The minimum is attained when (the two-photon amplitude (4) becomes a separable function in the idler and signal wave vectors). The Schmidt width for any of the above families of eigenmodes is always the same: . We stress that the Schmidt width does not represent the actual cross-section widths of the idler and signal beams. These widths, that can be measured experimentally, can also be obtained by resorting to the partially reduced density matrix of the idler and signal photons. For amplitude (4) one finds , which is consistent with the widths employed in (28) in the thin crystal approximation ().

Let us focus on the most usual encountered situation where the measured fiber modes are the fundamental Gaussian modes (of width at the phase plates). Since the involved eigenfunctions in the Schmidt decomposition (1) have cylindrical symmetry, we use the Laguerre-Gaussian modes as the convenient computational basis to derive the radial functions (see Appendix B). Remarkably, it turns out that the functions can be cast in terms of a sole parameter

(5)

where denotes the hypergeometric function, and

(6)

The parameter corrects by taking into account the spatial overlapping of the Schmidt and the fiber modes. Indeed, if and only if , does hold. When then , which corresponds to a constant phase-matching function. The functions increase monotonically from to for (). The fact that all depend on a single parameter will be exploited in the next Section to show how one can estimate the Schmidt number in a particular experimental scenario.

V Experimental Schmidt Number

Let us now examine the problem of measuring the amount of transverse entanglement of two-photon sources which, in our case, is characterised by the Schmidt number . In principle, a complete quantum tomography of the two-photon state could yield the desired amount of entanglement (11), but generally this would require a very large number of measurements, each for every possible pair of spatial modes, and this is technically demanding. Here we propose an alternative approach. In our simple model for the two-photon amplitude (4) we have found that the radial part of the coincidence probabilities (3) depends on a single parameter . This parameter involves the characteristic widths appearing in the two-photon amplitude (4), namely the pump beam width and the phase-matching width , together with the fiber mode width (at the phase plate locations). The phase-matching width could be measured by scanning in the plane of detection, but in our case it is not necessary. If, instead, one rotates the azimuthal phase plates (maintaining the detectors fixed), then, the recording of the coincidence distributions allows one to extract the value of as a fitting parameter for via Eqs. (3) and (5). Therefore, if is conceived as a parameter to be directly measured (rather than ), the amount of transverse entanglement can now be written in the form

(7)

where . This experimental Schmidt number depends on quantities that are easily accessible: the fitting parameter and the widths and (the presence of the width ratio should be interpreted as a correcting geometrical factor). It increases from , when , to infinity as .

By properly engineering the impulse response functions it is possible to enhance the sensitivity of with , thus improving the accuracy of the estimated . To this end, we consider two simple types of transparent azimuthal phase plates, and examine the dependence of Eq. (3) when the idler and signal phase plates are mutually rotated a relative angle . Their dispersionless impulse response functions are of the form . Owing to the initial perfect anticorrelation in OAM of the down-converted photon pairs, one expects that only when both photons are subjected to complementary azimuthal transformations the perfect anticorrelation is preserved and the coincidences are maximal. As soon as the phase plates are rotated in such a way that they are not longer oriented in a complementary arrangement, the photon coincidences decrease. Notice that the axes with respect of which the angles and are taken do not coincide. For symmetric phase plates these reference axes are inverted . In what follows, we assume that each plate is characterised by a noninteger parameter , where is the relative step height introduced by the plates, their refractive index, and the wavelength of the idler and signal beams.

Figure 2: (Left column) Normalised coincidence distributions as a function of the relative orientation between the idler and signal angular diaphragms. Solid, dashed and dotted curves are calculated for , and , respectively. (a) and (c) , both with . (Right column) The corresponding configurations of the angular diaphragms.

The first phase plate type we consider is an angular diaphragm [see inset in Fig. 2(a)]. It consists of a thin uniform dielectric circular slab with a “cake-slice” indentation that subtends an angle , with a nonzero (mod ) only if . Similar angular diaphragms have been employed in proof-of-principle demonstrations of the uncertainty relation for angular position and OAM (29). The probability (3) can be cast as . The main features of the normalized coincidence are: (i) it does not depend on the integer part of (it suffices to consider ), the visibility being maximal when ; (ii) the coincidence distributions are identical whether the aperture angle is or and symmetrical around ; (iii) for and fixed, the visibility diminishes as decreases; and (iv) the maximum visibility always occurs in the limit (constant phase-matching function and thus very high Schmidt number) where one has . Figures 2(a) and (c) depict the characteristic profiles of the (normalised) when the aperture angles of the angular diaphragms are and [Figs. 2(b) and (d) show, respectively, the configurations of the angular diaphragms that yield the coincidences (a) and (c)]. Comparing Figs. 2(a) and (c) one sees that the visibility of the former exhibits a much stronger variation with than the latter. This suggests that the phase plate configuration of Fig. 2(b) is preferred to that of Fig. 2(d) for achieving a more accurate estimation of the Schmidt number via Eq. (7). Finally, we should add that if is an integer number then the profiles of are constant (independent of ).

Figure 3: Normalised coincidence distributions as a function of the relative orientation between the idler and signal spiral phase plates. Solid, dashed and dotted curves are calculated for , and , respectively. (a) ; and (b) . (c) Configuration of the spiral phase plates.

The second type of azimuthal (refractive) component is a spiral phase plate [see inset in Fig. 3(a)] (19); (30). Its impulse response function includes the phase dependence , for . In this case, the joint probability (3) can be written as . At variance with , the coincidence distributions exhibit the development of interference ripples as decreases and increases [see Figs. 3(a) and (b)]. When a parabolic profile is obtained . In this limit the coincidences become independent of the integer part of , and, as with angular diaphragms, the maximum coincidence visibility is attained (when ).

We have considered other types of azimuthal phase plates and found the same parabolic dependence of the coincidences as . This is consistent with the abovementioned prediction that for very high transverse entanglement the joint probability becomes independent of the radial structure of the Schmidt modes. On the other hand, the absence of a vanishing minimum in the coincidence distributions is a signature of a finite amount of transverse entanglement. This effect gives rise to a contribution in the photocounts which always exists regardless of the detector efficiencies [see Fig. 2(a) and Figs. 3(a),(b)]. It can however be diminished for certain phase plate configurations [see Fig. 2(c)] and/or by employing a large ().

Figure 4: Coincidence distribution dependence on the relative orientation between idler and signal spiral phase plates (SPP with ). Circles are experimental values from Ref. (14), solid and dashed curves correspond to fitted with , and , respectively. The slight asymmetry in the experimental values was probably caused by the presence of a small anomaly in the center of the spiral phase plates (14). (Inset) Schmidt number as a function of the width ratio for .

Spiral phase plates were recently used in an elegant experiment aimed to show the high-dimensional spatial entanglement of a two-photon state from a down-conversion source (14). The relevant parameters of the experiment are: for the phase plates, m, pump width m, and thickness of the nonlinear crystal mm (14). According to the model presented in Ref. (31), which corresponds to our , it was concluded that . Figure 4 depicts the predicted distribution for (dashed line) together with the measured coincidences reported in Ref. (14). According to our approach, a probability distribution fitted with (solid line in Fig. 4) shows excellent agreement with the experimental results. Equation (7) then suggests that the corresponding Schmidt number should be much smaller than the above . However, we cannot perform a definite estimate for since the value of the fiber mode width at the location of the spiral phase plates, needed to calculate appearing in Eq. (7), was not provided in Ref. (14). The inset in Fig. 4 plots for a wide range of ratios . For instance, using our fitting parameter , together with m and assuming m, leads to and an estimate for , at least two orders of magnitude smaller than the quoted in Ref. (14). A very large value for could only be expected if, in the experiment, the fiber mode width in the plane containing the phase plates was much smaller than the pump width .

Vi Conclusions

In view of the situation described above, it is reasonable to conclude that additional experimental results and further theoretical analyses -a recent example can be found in Ref. (32)- are necessary and of interest in order to clarify the problem of quantifying the amount of transverse entanglement of photon pairs produced in down-conversion sources. If future values for the Schmidt number coming from accurate and solid measurements confirm the estimate of Ref. (14) and cannot be reproduced within our approach, one should conclude that the conventional Gaussian profile constitutes a poor approximation for the phase-matching function in Eq. (4). However, the ability of our predicted coincidence distributions to fit the experimental distributions of Ref. (14) suggests the correctness of our analysis and that the value for the Schmidt number quoted in Ref. (14) could be overestimated. In this case, our approach not only would provide a simple an accurate procedure to identify the relevant spatial modes and extract the degree of transverse entanglement; it would go in fact beyond by characterizing the action of all local bipartite azimuthal optical transformations on a broad family of two-photon states using only two detectors. Our results could also be of interest in other fields such as in identifying the intervening spatial modes when preparing well-controlled superpositions of photon states carrying OAM (to be coherently transferred onto Bose-Einstein condensates (33); (34) or by entangling photons with ensembles of cold atoms (35)), in high-resolution ghost diffraction experiments with thermal light (36), or in decoherence processes such as the entanglement sudden death mediated by the simultaneous action of several weak noise sources on bipartite photon systems (37).

Acknowledgements.
We thank S. Barnett, D. Diego, G. Molina-Terriza, M.J. Padgett and S. Walborn for useful discussions and gratefully acknowledge financial support from the Spanish Ministry of Science and Technology through Project FIS2005-01369, Juan de la Cierva Grant Program, CONSOLIDER2006-00019 Program and CIRIT Project SGR-00185.

Appendix A

In this Appendix we outline the derivation of the eigenvalues and eigenfunctions in the Schmidt decomposition (1) when the two-photon amplitude is given by Eq. (4). We start with some well-known facts.

For a general two-photon pure state, , it is possible to express as a bilinear sum of idler and signal basis states belonging to the Hilbert space of the system

(8)

where and label the set of quantum numbers for the idler and signal photons, respectively, whereas the coefficients describe the probability amplitudes for each tensor product of basis states. Our first aim is to evaluate the coefficients for the down-converted state . Let denote the basis wave functions in the transverse momentum representation. Inserting the closure relations , it can be readily seen that the coefficients are given by

(9)

Since we are interested in elucidating the correlation properties of the spatial degrees of freedom (OAM and the magnitude of the transverse radial wave vectors), we choose as the computational basis the complete set of normalised Laguerre-Gaussian modes (20)

(10)

where and denote the radial and azimuthal variables in momentum space, is the mode width (at the beam waist), and are the associated Laguerre polynomials. The indices and represent the winding (or topological charge) and the number of nonaxial radial nodes of the modes.

Combining Eqs. (4), (9) and (10), it follows that

(11)

By means of the well-known Anger-Jacobi identity , where is the modified Bessel function of the first kind, the two angular integrals yield the selection rule . This shows that the idler and signal photons are perfectly anticorrelated with respect to their topological charge, which is a manifestation of OAM entanglement. Hence, , where

(12)

The first of the radial integrals in (12) can be performed by resorting to the following formula

valid for all real , integers and , and complex ().

The second radial integral can also be done analytically. However, a dramatic simplification occurs if the widths and of the idler and signal radial modes, which at this stage have not been specified, are properly selected: . In this case, one obtains a second selection rule for the radial indices. Namely, . That is, with such a choice of the widths, which is unique, the idler and signal radial modes are perfectly correlated. But this shows that it is precisely for when one derives the Schmidt decomposition. Thus, the Schmidt width, , corresponds to . The coefficients reduce then to

(13)

Let . It is now clear from Eq. (13) that the Schmidt eigenvalues are . We therefore conclude that the Schmidt decomposition of the two-photon amplitude (4) is

(14)

where and represent the idler and signal Schmidt eigenmodes. These eigenmodes belong to the Laguerre-Gaussian basis. It is worth mentioning that although the Schmidt width is unique, other decompositions (actually infinitely many) are possible in different eigenmode bases. For instance, the Hermite-Gaussian modes (in Cartesian variables and the same ) constitute another possible basis. In fact, all eigenmode bases connected via the following unitary (metaplectic) transformation , where is a unit vector in the equatorial plane of the so-called orbital Poincaré sphere (parameterised by the spherical angles and ) and is an angular momentum operator, form Schmidt bases for (14). The components of the angular momentum operator: , , and are the three SU(2) generators. Of these, only represents real spatial rotations along the propagation of light, and it is thus the only component associated with the orbital angular momentum of light that can be measured in experiments (20).

Appendix B

In this second Appendix we explicitly show how to calculate the radial functions given in Eq. (5). For the chosen two-photon amplitude (4) and the measurement scenario, these radial functions depend on the spatial overlap of the Schmidt and the fiber (fundamental Gaussian) modes at the planes where the azimuthal phase plates are located. We need to evaluate the integrals . We use the Schmidt eigenvalues found in Appendix A. The idler and signal Schmidt modes and belong to the Laguerre-Gaussian mode basis. We stress once again the fact that the corresponding fiber and Schmidt widths are not necessarily equal. This implies that the orthogonality relation for these modes, , does not hold in general when .

We start by recalling a remarkable identity between Laguerre polynomials and the modified Bessel functions  (38)

(15)

valid for all real , integer and complex ().

In detail, the integrals to calculate read

(16)

Inserting identity (15) in (16) and employing the following result (38)

valid for and all complex , with representing the Whittaker functions, we obtain

(17)

where we have introduced the parameter . The remaining integral is carried out through the use of formula

valid when . The function denotes the hypergeometric function. Equation (17) is finally given by

(18)

This expression reduces to Eq. (5) when one leaves aside all the -independent factors (they do not play any significant role in the normalised coincidence profiles). In this way the two-photon detection probabilities (3) can be cast in terms of the impulse response functions , and the radial functions that solely depend on the parameter (when ).

Appendix C

In this appendix we would like to justify the model proposed here, the double gaussian profile as phase matching function.

a.1 Classical Energy for Non-linear Media

In this section we derive explicitly the energy of an anisotropic and non-linear medium. This classical description enables us to find a reasonable quantized Hamiltonian which allow us to describe the PDC process. In anisotropic and non-linear media we must be cautious in the description of physical variables. For a deep analysis, we begin by assuming that the polarization of the material obeys the usual expansion from non-linear optics

where (with ) represent the nonlinear susceptibility tensors of order responsible for the coupling of fields. The energy of the electromagnetic field can be written as

(20)

where the displacement vector field is defined as usual . In linear media, the energy density of the electromagnetic field depends quadratically on the electric field , whereas in non-linear (parametric) media, polarization exhibits a dependence on higher order electric field terms (LABEL:PDCPolarization). We can expand the energy (20) in two terms, one related to the linear contribution of the electric field and the other one related to the non-linear part. Taking the lowest order nonlinearity (LABEL:PDCPolarization), the energy due to this contribution can be written as (39); (40),

(21)

where the integration extends over the volume of the nonlinear medium. Notice that the second order susceptibility is not the same as the susceptibility in equation (LABEL:PDCPolarization). We will consider that the medium response to the electric field is not instantaneous. Hence, the second order polarization should be written as