Complex beam shaping by cascaded conical diffraction with intercalated polarization transforming elements.

Cascaded conical diffraction where optical elements modifying the local polarization state are intercalated between the aligned biaxial crystals is analyzed theoretically in the framework of paraxial diffraction theory. The obtained expressions are verified and confirmed experimentally for the case of a two-crystal cascade intercalated by a polarizer or a wave plate. The present approach can be used to realize a variety of vector beams with complex beam shapes composed of concentric rings with strongly modulated azimuthal intensity distribution. A potentially very fast switching of the overall beam shape is possible if the intercalated elements are electro-optically tunable retarders.


Introduction
When a light beam enters a biaxial birefringent crystal along one of its two optical axes it experiences a phenomenon known under the name of internal conical refraction (or internal conical diffraction). The beam propagates as a hollow cone inside the crystal and emerges as a diffracting hollow cylinder. When the proper observation plane is chosen (focal image plane) the transverse intensity distribution associated to this effect exhibits a double circular ring separated by a narrow dark region known as the Poggendorff dark ring. Even though this phenomenon was predicted by Hamilton [1] already in 1832 and was first observed by Lloyd [2] just one year later, investigations of conical diffraction experience presently a second life and a strong renewed interest both theoretically and experimentally , as recently reviewed by Turpin et al. [25]. This is due on the one hand to an improved understanding of the effect following its paraxial diffraction theory by Belskii and Khapalyuk [3] and its elegant reformulation by Berry in 2004 [4]. On the other hand the strong potential of this peculiar phenomenon for several modern photonics applications has been recognized. These include optical tweezers or bottle-type beams for trapping of particles [7][8][9], optical trapping of Bose-Einstein condensates [10], polarization metrology [11][12][13], polarization multiplexing for free-space optical communication [14] , super-resolution microscopy [15,16], lasers with specific polarization properties or spatial profiles [17][18][19][20][21], applications in the field of singular optics [22][23][24], and several other.
One of the major recent advances in the field of conical diffraction has been the extension of the effect to a cascade of two or more biaxial crystals with all their optical axes aligned. This approach adds versatility to the effect and leads in general to several concentric conical diffraction rings, with their relative intensities being governed by the angles of orientation of the crystals around the fixed direction of the optical axis. The general paraxial theory of cascaded conical diffraction was developed by Berry in 2010 [26]. An alternative approach based on the splitting and propagation of a bunch of classical rays was given by Turpin et al. [27,28] and several experimental and application oriented investigations of cascaded configurations were recently performed [8,14,23,24,[29][30][31][32][33]. In general the cascade of N crystals leads to 2 N −1 conical diffraction rings [26,28], for a circularly polarized or unpolarized input beam the intensity is azimuthally homogeneous on each of the rings. Notably, the local polarization on the rings is always linear with two radially opposite points exhibiting orthogonal polarizations, so that the conical diffraction process can be considered as a natural infinite channel polarization demultiplexer. These polarization properties suggest that the scrambling or filtering of the polarization between the crystals put in cascade should lead to a dramatic change of the overall observed conical diffraction pattern with respect to the case where the polarization is transferred without change. While the intercalation of polarization transforming elements in cascaded configurations was used in few experimental studies [8,24,29,32], a detailed theoretically description of this situation is still lacking.
In the present work we treat theoretically and experimentally cascaded conical diffraction where polarization transforming elements, such as wave-plates or polarizers, are intercalated between each pair of crystals. It is shown that the usual angular homogeneity of the intensity along the conical diffraction rings for unpolarized or circularly polarized input is lost. This leads to the possibility to realize complex vector-type light structures with highly localized distributions along the azimuthal direction. The involved intercalated optical elements are simple, spatially homogeneous and can be potentially switched very fast if realized with electro-optical devices. Therefore, intercalated cascaded conical diffraction can represent a valid and faster alternative to techniques based on pixellated phase elements (spatial light modulators, SLM) or liquid-crystal based q-plates to generate various classes of complex beam shapes [34,35]. Section 2 describes the terms of the problem and gives the theoretical treatment based on Berry's paraxial diffraction theory [26] extended to include the role of the intercalated elements. With the help of modified complex Belskii-Khapalyuk integrals we give explicit analytic solutions of the basic Fourier-type integral in the special case of a cascade of two crystals intercalated either by a λ/4-, a λ/2-plate or a polarizer. In Section 3 we give few specific examples and verify experimentally the theoretical predictions for the case of a cascade of two crystals of different length. Finally, the appendix gives some details on the method of solution of the basic Fourier integral that leads to the expressions given in Section 2 and defines the modified Belskii-Khapalyuk integrals.

Theory
We consider a series of N biaxial crystals arranged with a common optical axis oriented along the z direction so that each crystal, individually, gives rise to internal conical diffraction. As shown in Fig. 1, the crystals can be rotated with respect to each other around the common optical axis. The rotation of each crystal is expressed by the angle γ n between the x-axis and its direction γ n of displacement of the conical diffraction cone. Formally this direction is given by γ n ∝ k × ( k × S * ), where S * is the specific Poynting vector on the cone that gives the maximum walk-off angle with the wavevector k ||z (see right-hand inset in Fig. 1). Without loss of generality one can orient the first crystal parallel to the horizontal x-axis of the laboratory frame, so that γ 1 = 0. The orientation of the indicatrix (index ellipsoid) for this first crystal is indicated in the top left inset in Fig. 1. For the case γ 1 = 0 the longest and the shortest main axis of the indicatrix are in the laboratory xz-plane while the middle main axis is along the y-axis. For the following crystals the same projection of the indicatrix seen in Fig. 1 is found, however in a plane x z, where the axis x is rotated by an angle γ n around z with respect to the x-axis. Each crystal may be of a different length l n , also the crystals may be composed of different materials and thus be associated to different cone semi-angles α n . As done in [4,5,26] we consider a beam with an intensity 1/e radius w (measured in the focal plane of a lens placed before the crystals) and we normalize all transverse dimensions in real space with respect to this quantity w. In this way each crystal is characterized by a normalized strength parameter ρ n ≡ α n l n /w, which is the radius of the emerging ring in unit of the beam width. Between each pair of crystal we allow the presence of an optical element controlling the polarization state of the light, which is assumed to be oriented at an angle θ m with respect to the x-axis. These elements may be polarizers, quarteror half-wave plates, general wave retarders or combinations of any of these elements, they are characterized by a Jones matrix J m , where m extends to N − 1. Fig. 1. Arrangement for cascaded conical diffraction of N crystals with strength parameter ρ n intercalated by N − 1 polarization transforming elements with Jones matrices J m . The inset on the right shows the situation for the first crystal. The Poynting vector directions S associated to the common wave-vector k parallel to the optical axis lie on a cone containing the vector k. The direction γ 1 of displacement of the conical diffraction cone points towards the Poynting vector S * that has a maximum walk-off angle with k. The inset on top left shows the orientation of the projection of the index ellipsoid (indicatrix) on the xz-plane for the first crystal.
In order to fully account for diffraction we follow closely the Fourier optics approach reformulated by Berry [4] and include the effect of the polarization transforming elements. The optical beam entering the cascaded crystals is described as a sum of paraxial plane waves with wave-vector directions very close to the optical axis of the crystal. In Fourier space each wave-vector k = (k x , k y , k z ) can be represented in normalized cylindrical coordinates with the transverse components (κ, φ) such that k x ∝ κ cos φ and k y ∝ κ sin φ , where κ is normalized to 1/w, that is κ ≡ (k 2 x + k 2 y ) 1/2 w. Finally, we also use a normalized longitudinal spatial coordinate ζ [4,5,26] such that the position ζ = 0 corresponds to the plane where the rings are the sharpest. This is the "focal image plane" [4,5] of the incoming beam under the presence of all the crystals and all the polarization transforming elements. The normalization factor is given by the Rayleigh length z R = k 0 w 2 , that is ζ ↔ z/(k 0 w 2 ) with k 0 = 2π/λ and λ the light vacuum wavelength.
Let us place ourselves in the framework of the paraxial approximation and consider an input field for which the transverse distribution of the electric displacement vector is given in wavevector space by D 0 (κ, φ). The output distribution in real space D( ρ, ϕ, ζ ) can be obtained by Fourier transformation in polar transverse coordinates of this field, after being propagated through the whole optical system [26], that is In the common case where the input beam is homogeneously polarized and of circular symmetry the input field D 0 is expressed as where the function a(κ) gives the amplitude distribution as a function of the transverse wavevector and d 0 is a unit polarization vector. The matrix U tot in Eq. (1) gives the transfer function through the optical arrangement. In the case of a cascade of N conical diffraction crystals its expression was given in [26], where the U n (κ, φ, γ n ) are unitary matrices associated to the individual crystals and are expressed as In the case addressed in the present work, where the crystals are intercalated by polarization transforming elements, the matrix U tot should contain the effect of these elements, Eq. (3) should then be replaced by where the Jones matrices J m (θ m ) are not necessarily unitary. Obviously, in the absence of one or more of the polarization transforming elements the corresponding Jones matrices have to be replaced by the unit matrix. The intensity distribution is finally obtained up to an unimportant multiplicative constant from (1) and (5) as While the complex integral (1) can be determined by a rather lengthy brute-force numerical calculation, it is more convenient to perform the azimuthal integration analytically. Appendix A describes the method for the integration of the Fourier integral (1) and introduces the modified Belskii-Khapalyuk integrals B m ( ρ,ρ, ζ ). In the following we treat explicitly the specific cases where N = 2 and the intermediate element is either a quarter-wave plate, a half-wave plate or a polarizer. We limit ourselves to the case of a circularly polarized input wave, for which, in absence of intercalated elements, conical diffraction always leads to rings with azimuthally homogeneous intensity.

Two crystals intercalated by a λ/4-plate
We start by considering the special case where a λ/4-plate is placed between two crystals with parallel optical axes. The wave-plate is oriented under an angle θ with respect to the x-axis and the input polarization to the first crystal is homogeneous and circular so that In absence of the intermediate wave-plate this situation leads to two conical diffraction rings, each with a homogeneous intensity along the azimuthal coordinate. The intensity associated to each ring depends on the value of the orientation angle γ 2 of the second crystal with respect to the first [26,28]. The radii of the two rings in our normalized units are |ρ + | and |ρ − |, with The output light distribution in the presence of the λ/4-plate is calculated from (1) using the Jones matrix associated to the wave-plate. Using the modified Belskii-Khapalyuk integrals (21) defined in the Appendix A the two complex components D x and D y of the electric displacement vector and where we have used the abbreviation B m (ρ) ≡ B m ( ρ,ρ, ζ ). It can be easily seen from (8) and (9) that the resulting intensity distribution (6) is no longer independent from the real-space azimuthal angle ϕ, as will be discussed with the concrete examples in Sect. 3.

Two crystals intercalated by a λ/2-plate
In this case the relevant Jones matrix is and the procedure to calculate the output D-vector is similar as above. For homogeneously circularly polarized input one obtains instead of (8) and (9), and

Two crystals intercalated by a polarizer
Finally we consider the case where the intermediate element is a polarizer described by the non-unitary Jones matrix The components of the output D-vector are then and

Examples
In order to visualize the effects we give in this section some specific examples and compare them with corresponding experimental tests in the case of a two-crystal cascade. The experiments were performed at the wavelength of 633 nm using circularly polarized input light to the first crystal, which we chose to be the longest one. We used two crystals of KGd(WO 4 ) 2 (KGW) with lengths of 22.6 and 17.6 mm, respectively. By using the principal refractive indices n g , n m and n p of KGW determined by Pujol et al. [36], the cone aperture semi-angle α 1/(2n g n p )[(n 2 g − n 2 m )(n 2 m − n 2 p )] 1/2 is α 19.6 mrad for this crystal. For our focusing conditions with a f =100 mm spherical lens the beam 1/e half-width at the waist position is w 4.5 µm and the corresponding normalized strength parameters are ρ 1 98.6 and ρ 2 76.8.

Crossed crystals
Let us consider the case where the crystals are crossed with a relative angle γ 2 = π/2. In absence of any polarization transforming elements between them, this situation leads to two azimuthally  homogeneous rings (in fact two double rings), as shown in Fig. 2(a), obtained by integration of Eq. (1) for ζ = 0 in the case where the Jones matrix J 1 in (5) is identified to the unit matrix. As pointed out earlier [26,28], in this specific case the power is equally split among the two rings. Nevertheless the local intensity on the internal ring is larger due to the smaller area it occupies [28]. The latter is proportional to each ring radius and in our specific case the intensity ratio is roughly a factor of 8.
We consider first the intercalation of a quarter-wave plate oriented parallel to the first crystal (θ = 0). Figure 2(b) shows the experimentally observed intensity distribution in the two rings as obtained by imaging the focal image plane to a far away CCD camera by means of an imaging lens placed behind the crystal cascade. Figure 2(c) gives the corresponding expected intensity distribution calculated with Eqs. (8), (9) and (6) and with a(κ) = exp (−κ 2 /2) in Eq. (2) and the integrals (21). This distribution a(κ) corresponds to the Fourier spectrum of the input Gaussian beam associated to the 1/e half-width w. Figures 2(b) and 2(c) clearly show that the introduction of the λ/4-plate breaks the angular degeneracy and leads to an azimuthal dependence of the intensity in the two rings. In addition to the intensity distribution we plot in Fig. 2(d) also the distribution of the absolute value of the output displacement vector | D| (proportional to the square root of the intensity). Since this choice permits a better visualization of the weaker ring, we will keep this representation in the further examples. Figures 2(b), 2(c) and 2(d) clearly indicate that the intensity maxima and minima of the internal ring are in anti-phase with those of the external one. To look at this aspect in more detail we plot in Fig. 2(e) the expected intensities on the two rings as a function of the output angle ϕ. The blue solid line corresponds to the radius at which the internal ring has its maximum and the red dotted line is at the corresponding radius for the external ring. The corresponding experimental data are given in Fig. 2(f). The latter are obtained by evaluating the angular dependence of the average intensity within equally wide narrow rings that contain the internal and the external double rings, respectively. Figures 2(e) and 2(f) show a good agreement and confirm that the intercalated wave plate leads to two maxima and two minima for each ring, mutually in opposition of phase. The intensity modulation in the outer ring follows approximately a dependence in sin 2ϕ, while the one in the inner ring goes with − sin 2ϕ, this means that here the azimuthal intensity modulation is double as fast as in the case where linearly polarized light is used as input to the conical diffraction process [28]. In the case of a quarter wave plate the azimuthal intensity modulation contrast is approximately 1/2 for both the outer and the inner ring. With some little algebra one can show that this modulation contrast is roughly given by the ratio: 2 (Re[ 2 , where for each modified Belskii-Khapalyuk integral B m one should take B m ( ρ, ρ + , 0) for the outer ring and B m ( ρ, ρ − , 0) for the inner ring (see Appendix A). The above discussion holds for our example for which the λ/4-plate was oriented at the angle θ = 0. However, the choice of another angle θ leads solely to a rotation of the whole output intensity distribution by the double angle 2θ and all the conclusions remain therefore valid. We note also that, despite the use of the quarter wave plate, the light on the two rings is locally linearly polarized with a polarization direction depending on the angle ϕ. In the present case the output wave of the internal ring is horizontally polarized for ϕ = −90 deg, vertically polarized for ϕ = 90 deg and is polarized at +45 deg and -45 deg for ϕ = 0 and ϕ = 180 deg, respectively. Therefore, as is the case for conical diffraction in a single crystal, the local polarization direction angle is a linear function of ϕ/2 so that two opposite points on the same ring possess orthogonal polarizations. Also, for the same azimuthal angle the polarizations on the two rings are mutually orthogonal, which means that the polarization on the external ring for a given ϕ corresponds to the one on the internal ring for ϕ + 180 deg. We remark that the above polarization distribution is exactly the same as the one obtained for crossed crystals in absence of an intercalated element, i.e. the case of Fig. 2(a) for which the intensity is azimuthally homogeneous. Therefore, the intercalation of the quarter-wave plate does not modify the local polarization on the rings, this statement remains true also if the λ/4-plate is replaced by another polarization transforming element. The polarization distribution on the two rings depends only on the orientation of the two crystals.
It is worth noting that the intercalation of a half-wave plate according to Section 2.2 does not change dramatically the picture with respect to the above case of a quarter-wave plate and thus we discuss the differences only briefly. One gets also double peaked maxima and minima for each of the rings and, for a same orientation of the wave plates the intensity distribution keeps the same overall orientation. However, for the case of a half-wave plate the modulation is complete for both the inner and the outer ring. For example, if θ = 0 one gets zero intensity points for the internal ring at ϕ = π/4 and 5π/4 and zero intensity points for the external ring at ϕ = −π/4 and 3π/4. Therefore the use of a variable retardation wave plate permits to tune the azimuthal intensity modulation from zero to full contrast going from zero retardation to half-wave retardation, and back to zero if the retardation ranges between λ/2 and a full wave.
As a next example we consider the case where a polarizer under the angle θ = π/2 is inserted between the crossed crystals. As shown in Fig. 3 this leads to a richer and more complex azimuthal intensity redistribution within each of the two rings. Figure 3(a) shows the expected distribution of the modulus of the D vector as obtained from Eqs. (14) and (15). It clearly shows that the intensities in each of the two rings exhibit a well identified maximum at a certain angle and that the mutual positions of these maxima are angularly shifted by ∆ϕ = π/2. This is clearly seen also in the experimental observation of Fig. 3(b) and in the theoretical and experimental intensity distribution along the two rings depicted in Fig. 3(c). Each ring possesses two zero-intensity points with a weak secondary intensity maximum between them, as seen for instance in the inset in Fig. 3(c). The zero-intensity points form a right angle with the center of the cone projection, in our specific case they are found at ϕ = −π/2 and ϕ = 0 for the internal ring, and at ϕ = 0 and ϕ = +π/2 for the external one. The 90 degrees out-of-phase relative orientation of the internal and external rings may be interpreted as the manifestation of some kind of "pseudo-chirality" associated with the two-crystal structure. In the case discussed here, for which γ 2 = +π/2, this chirality is positive as defined by the fact that the external ring is oriented 90 degrees clockwise with respect to the internal one. The reverse would be true if the direction of the second crystal is inverted and γ 2 = −π/2. Interestingly, such an inversion of chirality is obtained also if the order of the crystals is reversed and the shorter crystal would be put before the longer one. This is in contrast to the above case where the intercalated element is a wave plate, for which the overall orientation of the rings is independent of the order in which the birefringent crystals are put into the set-up. This difference is associated to the fact that in the case of a polarizer the total energy in the beam is not conserved and the corresponding Jones matrix (13) is not unitary. Finally, as was the case for the intercalated wave plates, it is worth noting that rotating clockwise the polarizer by an additional angle ∆θ leads to a clockwise rotation of the whole light distribution structure by the double angle 2∆θ. For instance, the structure for a polarizer under the angle θ = 0 is obtained by central point symmetry from the one depicted in Fig. 3.

Parallel crystals
We consider now the interesting situation where the two crystals are parallel with the vectors γ 1 and γ 2 oriented in the same direction, so that the relative angle is γ 2 = 0. In this case, in absence of intercalated elements one obtains a single double ring at ρ ≈ ρ 1 + ρ 2 [26,28]. This radius corresponds to the one of the external ring in Fig. 2(a). The system behaves like there was only a single crystal having a length corresponding to the sum of the lengths of the individual crystals.
We first consider the case where the intermediate local polarization is scrambled by a quarterwave plate. As seen in Fig. 4 clearly the presence of the wave plate "reactivates" the otherwise non existing internal ring. As for the case of the crossed crystals given in Fig. 2 the intensities are azimuthally modulated along each of the rings with a two-fold symmetry. The maxima and minima of each ring are again mutually out-of-phase. However, here the modulation depth of the internal ring is complete, while the one of the external ring is only partial like in the case of the crossed crystals. Again, as in the previous cases, any further rotation ∆θ of the wave plate leads to a rotation of the whole light distribution structure by 2∆θ. It should also be noted that the use of a half-wave plate instead of a quarter-wave plate leads to a qualitatively similar and equally oriented picture, the only major change being that not only the internal, but also the external ring has a full modulation contrast. The fact that the internal ring can be activated by the presence of the wave plate has potentially very interesting applications. For instance, the use of an electro-optically tunable retarder shall allow to switch on and off very rapidly the internal ring structure without the need for any moving parts. The same is true for the switching on and off of the external ring in the case where the two crystals are oriented in antiparallel direction (γ 2 = π). As a final experimental example we discuss briefly the case where a polarizer is inserted between the parallel crystals. As seen in Fig. 5 also here the internal ring is reactivated with the same two-fold symmetry as for the case of Fig. 4. However, in contrast to all previous cases, the overall rotational symmetry associated to the intensity in the external ring differs from the one of the internal one and is only one-fold.

More than two crystals
Finally we give briefly two examples for the more complex cases where three or four crystals are put in cascade in the way shown in Fig. 1. In principle also for N > 2 a semi-analytical treatment like the one in Sections 2.1 to 2.3 can be done, however the related expressions become quite lengthy. It can be easily shown that for a number N of crystals put in cascade and intercalated by various polarization transforming elements the maximum order for the integrals B m involved in such expressions is m = N. Therefore for three crystals the integrals B 3 ( ρ,ρ, ζ ) are needed in addition to B 0 , B 1 and B 2 , while for four crystals also the B 4 ( ρ,ρ, ζ ) are required. Here, instead of using the semi-analytical approach, we show in Fig. 6 the expected intensity distributions obtained by a direct numerical integration of Eq. (1) in the focal image plane ζ = 0 and for a circularly polarized input wave.  Fig. 6 is for the case of three cascaded crystals with the first two crossed to each other and the last two parallel to each other. A quarter-wave plate oriented at 45 deg is placed between the first two crystals and a half-wave plate oriented at θ 2 = 90 is placed between the last two crystals. The chosen normalized strength parameters of ρ 1 = 40, ρ 2 = 15, ρ 3 = 5 lead to four conical diffraction rings for which the normalized radii are roughly ρ ≈ (20; 30; 50; 60), which can be easily recognized in Fig. 6(a). The azimuthal intensity distribution on each of the rings exhibits two nodes which are aligned horizontally for the most internal (ring 1) and the most external ring (ring 4), and vertically for the two intermediate rings. Remarkably, in the present configuration the azimuthal intensity profile between these nodes is not symmetric, so that the intensity center of mass shifts clockwise for the rings 1 and 3, and counterclockwise for the rings 2 and 4.

Panel (a) in
For the case of Fig. 6(b) a fourth crystal with ρ 4 = 70 and γ 4 = 135 deg is added to the previous three. Here we consider the case where polarizers are inserted between the outer crystal pairs (θ 1 = 0 and θ 3 = 135 deg) and a half-wave plate is put between the second and third crystal (θ 2 = 45 deg). The expected radial positions of the eight resulting conical diffraction rings are ρ ≈ (10; 20; 40; 50; 90; 100; 120; 130). All these rings can be recognized in Fig. 6(b), however the ring at ρ ≈ 90 is only hardly visible due to a very weak associated intensity. As was the case for instance in Fig. 3, the polarizers lead to a quite complex angular dependence of the intensity on each of the eight rings and the angle ϕ associated to the maximum intensity differs for each of them. Also, the increase of the number of crystals and intercalated elements leads to a sharper confinement of the intensity on a rather narrow angular region for the main lobe of each ring.

Conclusions
We have described theoretically and verified experimentally the effect of intercalating polarization transforming optical elements between the biaxial crystals forming cascaded conical diffraction. The additional elements break the usual azimuthal intensity homogeneity expected for unpolarized or circular polarized input beams. Complex vector-type beams are obtained with their shapes governed by the crystals conical diffraction strength parameters ρ n , their angular orientations γ n and the nature and orientation of the polarization transforming elements. A particularly interesting case is the one where otherwise silent rings are "re-activated" by the presence of the intercalated elements, which occurs if two or more crystals are arranged parallel or anti-parallel to each other. Since variable retarders can be realized by means of electro-optical devices, this opens the possibility to switch on and off individual conical diffraction rings at speeds exceeding several MHz. We have given explicit analytic expressions only for the case of a two-crystal cascade with intercalation of a λ/4-plate, a λ/2-plate or a polarizer. However, the general formalism remains valid also for other polarization transformations between the crystals, as may be obtained for instance by a combination of optical elements described by an appropriate Jones matrix. Our few examples have shown that a variety of complex beam shapes with strong (and different) azimuthal localization of the light intensity on each ring can be obtained. Obviously this localization could be improved even further by post-filtering the polarization state after the last crystal by means of a polarizer. The richness and versatility of the vector beam shaping features resulting from the present approach open up interesting perspectives for virtually every application that has been proposed in connection with conical diffraction and complex beam shaping, including fast switchable optical trapping, singular optics, material processing, polarization metrology, and super-resolution microscopy.

Appendix A: Integration of the Fourier integral (1)
When evaluating the form of the 2×2 complex matrix U tot (5) or of the product U tot · d 0 inside the integral (1) one gets a complicated sum over sine and cosine functions containing as arguments various linear combinations of the angles φ, θ n , γ n as well as the products κ ρ n . It is convenient to express this trigonometric sum as a sum of exponential functions of the same arguments. The general form of an individual term may then be expressed as where m is a positive integer and β is a linear combination of the angles θ n and γ n . The quantitỹ ρ is composed of simple sum or differences of the normalized strength parameters ρ n . For instance in the case where we have only two crystals (N = 2) the possible values ofρ are ρ 1 + ρ 2 , ρ 1 − ρ 2 , ρ 2 − ρ 1 and − ρ 1 − ρ 2 . The values ofρ are associated to a specific conical diffraction ring. However, sinceρ and −ρ belong to the same ring, the total number of concentric rings observed in cascaded conical refraction is only 2 N −1 and not 2 N , as pointed out earlier [26,28]. Considering only the azimuthal φ-integral in (1) for the specific term (16) we get The above integral is a special form of an integral of the following class for which a general solution was given by Massidda [37], 1 2π 2π 0 e imφ e a cos(φ −α 1 ) e 2b cos 2 (φ −α 2 ) dφ = e b e imα 1 k e 2ik (α 1 −α 2 ) I 2k+m (a)I k (b) , (18) where I k (x) is the modified Bessel function of the first kind of order k, which is related to the Bessel function of the first kind J k (x) by I k (x) = (1/i k ) J k (ix). Since in (17) b = 0, only the term k = 0 contributes to the sum on the right-hand side of (18). With a = iκ ρ, α 1 = ϕ and using Therefore, upon integration the azimuthal phase e ±imφ in Fourier space leads to a corresponding phase e ±imϕ containing the azimuthal angle ϕ in real space. One can now use the integral p to evaluate the contribution q of the term (16) to the integral (1), one obtains, The quantities B m ( ρ,ρ, ζ ) in the above expression are modified Belskii-Khapalyuk integrals that we define as Note that, unlike for the standard Belskii-Khapalyuk integrals [5,38], the above modified integrals are complex even in the focal image plane ζ = 0. However, it follows directly from the above definition that for this plane the following symmetries hold Fig. 7 visualizes the real and imaginary parts of the modified Belskii-Khapalyuk integrals B 0 , B 1 and B 2 for the case ζ = 0. It can be easily recognized that in this case the integrals assume significant values only in the neighborhood of the normalized radius ρ = |ρ|, which is indicated by the vertical lines in Fig. 7. The radii ρ = |ρ| withρ = ρ 1 + ρ 2 orρ = ρ 1 − ρ 2 correspond roughly to the radii of the Poggendorff dark rings of the two-crystal cascaded conical diffraction. Note that the rather sharp curves for B 0 , B 1 and B 2 in Fig. 7 are due to the choice ζ = 0. If we leave the focal image plane (ζ 0) the curves become much broader and can assume significant values even for ρ being far from the characteristic valuesρ, reflecting the fact that the conical diffraction rings get defocused. Note also that in the specific example shown here, which is associated to large values ofρ, the curves for B 0 , B 1 and B 2 appear very similar, even though they are not identical. The differences between the curves become much more pronounced for smaller values ofρ (not shown in Fig. 7), as obtained for a less focused input beam, for shorter crystals or for materials with a smaller aperture angle of the conical diffraction.
Finally we remark that in the analytic expression for the output components of the electric displacement vector given in sections 2.1 to 2.3 one always finds sums (or differences) of the modified Belskii-Khapalyuk integrals for values ofρ of opposite sign. Specifically these are B 0 (ρ ± ) + B 0 (−ρ ± ), B 1 (ρ ± ) − B 1 (−ρ ± ) and B 2 (ρ ± ) + B 2 (−ρ ± ) (see Eqs. (8) and (9) as well as the corresponding equations in sections 2.2 and 2.3). In combination with the symmetries expressed by Eqs. (22) and (23) this implies that only the real parts of the integrals B 0 , B 1 and B 2 play a role for the intensity distribution in the plane ζ = 0, what is no longer true if one leaves this plane. Also, even though we prefer to stick to the general form of the modified integrals (21), it is worth mentioning that the above sum (or differences) are directly proportional to the standard form of the Belskii-Khapalyuk integrals. In the latter, instead of the complex term exp(iκρ), the integrand contains a term cos(κρ) for m = 0 and for its generalization to all even values of m, and contains a term sin(κρ) for m = 1 and for the generalization to all odd values of m.  (21) as a function of the normalized radius ρ for the case ζ = 0 and a(κ) = exp(−κ 2 /2). The panels in the left column give the real part and the panels in the right column give the imaginary part of the integrals. The top panels are for B 0 (ρ + ) (solid lines) and B 0 (ρ − ) (dotted lines). The corresponding functions for B 1 (ρ ± ) and B 2 (ρ ± ) are in the middle panels and bottom panels, respectively. Hereρ + ≡ ρ 1 + ρ 2 andρ − ≡ ρ 1 − ρ 2 , with ρ 1 = 98.6 and ρ 2 = 76.8. The vertical lines correspond to the conditions ρ =ρ − and ρ =ρ + .