## Abstract

Using a system of coupled nonlinear Schrödinger equations (CNLSEs), we show that nonlinear light propagation in self-focusing Kerr media can be controlled via a suitable combination of linear and circular birefringences. In particular, magneto-optical effects are taken as a specific physical example, which enables the introduction of both types of birefringences simultaneously via the joint action of the Cotton-Mouton and the Faraday effect. We demonstrate the efficient management of the collapse of (2 + 1)D beams in magneto-optic dielectric media, which may result in either the acceleration or the suppression of the collapse. However, our study also shows that a complete stabilization of the bimodal beams (i.e., the propagation of two-dimensional solitary waves) is not possible under the proposed conditions. The analysis is performed by directly numerically solving the CNLSEs, as well as by using the variational approximation, both showing consistent results. The investigated method allows high-power beam propagation in Kerr media while avoiding collapse, thus offering a viable alternative to the techniques applied in non-instantaneous and/or non-local nonlinear media.

© 2010 OSA

## 1. Introduction

The nonlinear Schrödinger equation (NLSE) is the ubiquitous model describing wave propagation in nonlinear media. It applies to various physical contexts in optics, hydrodynamics, condensed matter and plasma physics, as well as in other fields [1–6]. In different settings (in particular, for different spatial geometries), the solutions of the NLSE depict different phenomena, such as the generation of stable solitons and self-similar structures, as well as wave collapse. Particularly, in the (1 + 1)D case, when only one transverse spatial dimension is involved in the dynamics, the NLSE is integrable by means of the inverse scattering transform, giving rise to stable solutions that preserve their form after interactions and collisions (i.e., solitons) [7]. Higher-dimensional NLSEs are not integrable, and hence no exact analytical solutions are known. Specifically in optics, stationary solutions to the (2 + 1)D NLSE with a self-focusing nonlinearity in two transverse dimensions do exist in the form of axially-symmetric *Townes solutions* [8–10], which are restricted to a single (*critical*) value of the beam power. These solutions are unstable, eventually leading to either a catastrophic collapse (singularity) or diffraction. Moreover, when the power exceeds the critical value, wave collapse always occurs [3,10,11], resulting in significant changes in the spatial profile of the beam (e.g., its filamentation), or alternatively, in material damage in real physical situations. Another example of such a catastrophic self-focusing is found in self-attractive Bose-Einstein condensates (BECs), when the total number of atoms is larger than the critical value [12,13].

Following its prediction in the 1960s [14], many experimental and theoretical works have been devoted to further investigations of the collapse dynamics in various systems (see e.g., [9–45]), including bulk self-focusing Kerr media [15,16], BECs [12,13], fluids [17] and plasmas [15,18]. An important subject in this field of interest is the possibility to *control* and eventually *arrest* the collapse. In this context, the general case of the stabilization of spatial-temporal collapse, yielding to the formation of light bullets (i.e., non-dispersing, non-diffracting optical pulses) was predicted in [21] and intensively investigated (e.g., see [22–24] and references therein). Various physical mechanisms to delay or even to completely eliminate the collapse have been suggested and described theoretically by properly modifying the NLSE [23,25]. In particular, it has been shown that the collapse can be arrested by including additional physical effects, such as multiphoton absorption and plasma formation [26], temporal dispersion (associated to the propagation of ultrashort pulse) [25,28], energy dissipation [25], changes in the propagation constant throughout propagation [29], and nonparaxial corrections [25,27,30]. Another mechanism, which allows direct arrest of the collapse, can be provided by the discreteness of the medium [31] and for instance it has been implemented by the replacement of a bulk medium by a bundle of parallel fiber-like waveguides, such as those that may be written in silica by means of femtosecond laser ablation [32]. On the other hand, the partial spatial coherence of the beams propagating in inertial bulk Kerr media has been shown to increase the critical power, thus delaying the collapse, but without the possibility of arresting it completely [33]. Recently, it has been demonstrated that the spiraling of self-trapped beams is an efficient mechanism for suppressing collapse in pure Kerr media (supporting stable elliptic solutions) [34]. Other approaches allowing for the stabilization of the catastrophic self-focusing rely on the competition between a quadratic (χ^{(2)}) and a cubic (χ^{(3)}) nonlinearity [35], as well as on the use of saturable [25,36,37] and nonlocal [25,38,39] nonlinearities. In optical systems, the nonlocality may be achieved in photorefractive materials or in liquid crystals, allowing for the formation of stable multidimensional spatial solitons [e.g., 38,39]. Specific types of nonlocal attractive interactions in BECs (in particular, dipole-dipole interactions) have also been shown to be suitable in arresting the collapse and stabilizing (2 + 1)D solitons [40]. Moreover, it has been demonstrated that (2 + 1)D quasi-stationary solitons may be stabilized if the beam is launched through a suitable periodic structure consisting of alternating layers of self-focusing and defocusing materials [41]. A counterpart of the latter method has been suggested in BECs, based on the use of the Feshbach-resonance technique in a variable external field [42,43]. Recently, passive optical elements has been proposed as another possibility for controlling the collapse of laser pulses in air [44]. On the other hand, both fundamental solitons and solitons carrying an embedded vorticity can be readily stabilized by transverse lattices, inducing a spatially-periodic potential that may suppress both their collapse and their decay [28,45].

In this work, we study theoretically a possible method for the management of collapsing bimodal optical beams by means of specially adjusting the birefringences, which in turn allow one to control the energy and phase transfer between the two-orthogonal-polarization components of the beam. The nature and the fundamental character of the interplay between birefringence and nonlinearity have been previously demonstrated for temporal solitons in optical fibers [46–49]. Here, using both the direct simulations for coupled (2 + 1)D NLSEs, as well as the variational approximation (VA) applied to the same system, we demonstrate that optical birefringence can be used as a powerful tool for the control of the onset of the collapse of two-component beams in bulk self-focusing Kerr media. Specifically, a system of (2 + 1)D CNLSEs has been previously applied to determinate the effect of four-wave-mixing (FWM) and walk-off [50,51], as well as to analyze the influence of the degree of coherence [33] on the critical power for which collapse takes place in Kerr media. It has been demonstrated that FWM contributes to self-focusing by lowering the threshold for collapse, whereas the walk-off and partial coherence act against collapse (i.e., delay it) by increasing the critical power. In the considered case, we propose to take advantage of the combined effects of the linear and the circular optical birefringences, which can lead to the acceleration, the delay, or even the effective suppression of the optical collapse. In particular, we analyze the interplay [52] of the Cotton-Mouton (CM) [53,54] and Faraday [53,55] effects, in order to generate the necessary combination of the magnetically-induced birefringences [56] that may control the collapse caused by the self-focusing Kerr nonlinearity in magneto-optic (MO) bulk media. However, because birefringence is a linear effect, which simply allows an energy exchange between beam polarization components, it does not lead to the full stabilization of high-power beams, when acting alone (as we demonstrate here). Nevertheless, the proposed technique is capable to significantly delay the onset of collapse, similarly to the birefringence-induced delay of the supercontinuum generation [57].

## 2. The theoretical model

#### 2.1. Coupled nonlinear Schrödinger equations (CNLSEs)

Our theoretical analysis is based on a system of coupled NLSEs describing the evolution of the scaled complex electric-field amplitudes *u _{L}* and

*u*of the left- (

_{R}*L*) and right- (

*R*) circularly polarized components, in the presence of the Kerr nonlinearity [5] and under the influence of the combined (here magnetically-induced) linear and circular birefringences [46–49]. By assuming an axisymmetric configuration, and limiting the analysis to states with zero vorticity (i.e., fundamental beams), one may assume that the solutions of the propagation equations take the form of ${u}_{L,R}\left(x,y,z\right)={u}_{L,R}\left(r,z\right)$, where $r\equiv \sqrt{{x}^{2}+{y}^{2}}$is the radial coordinate in the transverse plane. Assuming the validity of the paraxial approximation, the respective partial differential equations (PDEs), written in terms of rescaled coordinates $\tilde{z}$ and $\tilde{r}$, take the following form:

Typical normalization factors were used to derive the equations in their dimensionless form. Namely, we set $\tilde{r}=r/{w}_{0}$, $2\text{\hspace{0.17em}}\tilde{z}=z/{L}_{D}$ (note that the tildes are suppressed in the derivations of the variational approximation presented below) and ${u}_{L,R}={k}_{0}{w}_{0}{n}_{0}{U}_{L,R}\sqrt{{n}_{2}/(3\eta )}$. Here *w*_{0} is the initial radius (waist) of the axisymmetric Gaussian beam, ${L}_{D}=\pi {w}_{0}^{2}{n}_{0}/{\lambda}_{0}$ is the diffraction length, and *U _{L}*

_{,}

*are the complex amplitudes of the two components of the monochromatic field [i.e., the full components corresponding to the circular polarizations are defined as ${E}_{L,R}\left(r,z,t\right)={U}_{L,R}\left(r,z\right)\text{\hspace{0.17em}}\mathrm{exp}\left(i{n}_{0}{k}_{0}z-i\omega t\right)$]. In addition,*

_{R}*k*

_{0}is the wave number in vacuum, λ

_{0}is the wavelength,

*n*

_{0}is the linear refractive index,

*n*

_{2}is the nonlinear Kerr coefficient, and η is the vacuum impedance. Equations (1) were derived by means of the slowly-varying envelope approximation (SVEA), and also neglecting the propagation losses. Indeed, in the materials typically used in the experimental conditions, such as YIG crystals, the latter postulation is acceptable, due to the intrinsic properties of the medium [61]. We also assumed that the nonlinear and diffraction length-scales are much shorter than the dispersion length, thus allowing us to neglect temporal-dispersion effects.

The ratio of 2 between the coefficients in front of the XPM and SPM (cross- and self-phase modulation) terms in the nonlinear part of Eqs. (1) is a well-known feature of the circular-polarization basis, in the case of cubically-symmetric Kerr media. This basis has been chosen as the most convenient for the analysis of the considered problem. An equivalent physical model represented via the linear polarization states may also be analyzed, but this is not necessary, as its solutions are equivalent to those derived from Eqs. (1). Note that when the equations are written for the linear polarizations, four-wave mixing terms explicitly appear (see e.g., [50,51]) rendering the numerical analysis more complicated.

Of particular importance are the terms in Eqs. (1) with real coefficients *b* and *c*, which (in the considered case of MO materials) depend on the strength and the direction of an external magnetic field. By definition, both coefficients *b* and *c* are positive, and they represent the magnetically-induced circular and linear birefringences, respectively. Namely, *b* accounts for the Faraday effect, while *c* accounts for the Cotton-Mouton effect. These scaled coefficients are related to their physical counterparts, $\Delta {n}_{F}$ and $\Delta {n}_{CM}$, as follows: $b={k}_{0}^{2}{w}_{0}^{2}{n}_{0}\Delta {n}_{F}/2$ and $c={k}_{0}^{2}{w}_{0}^{2}{n}_{0}\Delta {n}_{CM}/2$. It is evident from Eqs. (1) that the linear and circular birefringences determine the rates of the linear amplitude mixing and phase shift between the two circularly-polarized fields, respectively.

To conclude the description of the model, it is worth mentioning that by replacing the propagation coordinate *z* with the time coordinate *t* in Eqs. (1), one can obtain the set of coupled Gross-Pitaevskii equations (GPEs) that is commonly used to describe BEC mixtures of two different spin states of the same atomic species [58].

#### 2.2. The variational approximation (VA)

In addition to direct numerical solutions of the CNLSEs, the variational approximation is also suitable for the description of the beam propagation and significantly reduces computation time. In our case, the VA amounts to a system of ordinary differential equations (ODEs) for the evolution of the width, amplitude, and relative phase of the two polarization components [59]. It will be demonstrated in the next section that the VA provides essential and accurate information about the beam evolution throughout propagation. In the present model, the VA yields reasonable predictions about the collapse dynamics, and indicates a specific range of parameters (input optical powers and birefringence coefficients) where more detailed, time-consuming direct numerical simulations of the CNLSEs are preferable.

Equations (1) can be derived from the Lagrangian $L=2\pi {\displaystyle \int \Lambda \text{\hspace{0.17em}}r\text{\hspace{0.17em}}dr}$ with density Λ:

Assuming, as above, an axially symmetric beam, the solutions of Eqs. (1) may be approximated by the Gaussian ansatz:

*A*,

*φ*,

*θ*,

*w*

_{R}_{,}

*and*

_{L}*ξ*

_{R}_{,}

*are real functions of*

_{L}*z*. Specifically,

*A*and

*φ*are the common amplitude and phase of both polarization components,

*θ*determines the distribution of the power between them,

*ψ*is the relative phase, and

*w*

_{R}_{,}

*and*

_{L}*ξ*

_{R}_{,}

*are the widths and radial chirps, respectively.*

_{L}The general ansatz in the form of Eq. (3) leads to cumbersome variational equations. A more tractable form of the approximation is based on a simplified symmetric trial function with ${w}_{R}={w}_{L}\equiv w$ and ${\xi}_{R}={\xi}_{L}\equiv \xi $, i.e., assuming the same widths and chirps for both components:

The total scaled power (conserved during propagation) corresponding to ansatz (4) is:

The substitution of Eq. (4) in Eq. (2) and the integration in the transverse plane yields the effective Lagrangian:

The Euler-Lagrange equations, in the general form of $\delta {L}_{eff}/\delta \zeta =\partial {L}_{eff}/\partial \zeta -\partial (\partial {L}_{eff}/\partial \zeta )/\partial z=0$ (where *ζ* is any variational parameter), are used to derive the evolution equations from the Lagrangian, Eq. (6). The first equation, $\delta {L}_{eff}/\delta \varphi =0$, reduces to the trivial solution $d\tilde{N}/dz=0$, which implies conservation of the total power. The second equation, $\delta {L}_{eff}/\delta \xi =0$, yields a well-known expression for the chirp, $\xi ={w}^{\prime}/w$. This relation may be combined with the equation $\delta {L}_{eff}/\delta w=0$ to eliminate *ξ* and to derive a second-order ODE for the evolution of the beam width *w*(*z*):

The remaining equations, $\delta {L}_{eff}/\delta \psi =\delta {L}_{eff}/\delta \theta =0$, give rise to the following subsystem of ODEs:

Equations (7) form the final set of the dynamical equations describing the evolution of the beam in the framework of the VA. The corresponding Hamiltonian (which is another dynamical invariant) is:

Incidentally, we note that by making use of the Hamiltonian of the system, it may be interesting to try to derive the analytical criteria for the collapse via an appropriate generalization of the virial theorem [15]. The latter may be useful to estimate the threshold power and identify regions of different dynamical behaviours associated to beam propagation (e.g., [33,34,50,51]). However, this analytical approach can be found to be inaccurate (i.e., in general determined conditions should be verified numerically) and as such, we have decided not to use it in the present work.

When analyzing the dynamical equations it is worth to note that in the case of *c* = 0 (with no linear mixing between the circular polarization components), Eq. (7b) yields *θ* equal to a constant. Concurrently, Eq. (7a), which governs the evolution of the beam width, decouples from the other equations (as expected; see [60]). In particular, for *θ* = 0 it can be demonstrated, using Eq. (7a), that, when $\tilde{N}<{\tilde{N}}_{cr}$, all solutions behave in such a way that $w(z)\to \infty $ at $z\to \infty $, thus implying the decay of the beam. On the other hand, the collapse (for which $w(z)\to 0$ at a finite value of *z*) always takes place when $\tilde{N}$ exceeds the critical value${\tilde{N}}_{cr}=1$. In fact, the latter condition is a well-known variational prediction for the critical power of *Townes* solitons [60]. For $c\ne 0$, the evolution of *w*(*z*), as governed by
Eq. (7a), no longer decouples from the evolution of *θ*. In fact, while the term sin^{2}(2*θ*) may take values in the range between 0 and 1, ${\tilde{N}}_{cr}$formally varies within the interval defined as $2/3\le {\tilde{N}}_{cr}\le 1$, thus suggesting a possibility for the stabilization against the collapse.

In addition, two other particular ranges of parameters can be considered for the analysis of the dynamical evolution equations. In the limits of large *c*, Eq. (7a) takes the form of a non-autonomous ODE decoupled from Eqs. (7b) and (7c):

*b*, is large, Eq. (7a) simplifies into the following form:

It is important to note that Eq. (10) includes a small coefficient in front of the cosine term, which is not present in Eq. (9), hence it is expected to lead to a much smaller stability region than in the case of large *c*.

#### 2.3. Magnetically-induced birefringences

MO effects, i.e., the modification of the optical properties of dielectric media by an external magnetic field, allow one to induce a certain birefringence in a medium that is otherwise optically isotropic (in the absence of an external magnetic field). This can be described by a proper modification of the dielectric permittivity tensor, *ε*, whose elements depend on the external magnetic field. Assuming that a magnetic field ** H** is applied in the (

*x*,

*z*) plane at an angle

*α*with respect to the

*z*-axis, the dielectric permittivity tensor is written as [52,53]:

*n*

_{0}is the refractive index in the absence of the magnetic field, and

*Q*,

*B*

_{1},

*B*

_{2}are material parameters depending on the magnetization of the medium, which vanish in the absence of the field (in particular,

*Q*~|

**| and**

*H**B*~|

_{i}**|**

*H*^{2}). The components of the dielectric tensor which are quadratic in |

**| account for the CM effect, while those linear in |**

*H***| represent the Faraday effect in transmission, or the MO Kerr effect in reflection (the latter one should not be confused with the commonly-known Kerr effect which gives rise to the electric nonlinearity of dielectric materials). In general, if the diagonal terms of the permittivity tensor are different, the material is linearly birefringent, whereas the off-diagonal terms, ~**

*H**iQ*, imply the presence of an optical activity (circular birefringence). Assuming that the beam propagates along the

*z*-axis, and defining the elements of the dielectric tensor as: ${\epsilon}_{1}={n}_{0}^{2}+{A}_{1}{\left|H\right|}^{2},$$\text{}{\epsilon}_{2}={n}_{0}^{2}+{A}_{2}{\left|H\right|}^{2},$${\epsilon}_{3}={A}_{3}\left|H\right|$, where

*A*are appropriate material constants, the following equation describes the resulting refractive index

_{i}*n*(square root of the dielectric permittivity) of the medium under the influence of the total magnetic field

**[52]:**

*H*The general solution of Eq. (12) is:

For each value of the angle *α* two different solutions for *n*^{2} exist, indicating that the intrinsically isotropic medium becomes *birefringent* under the action of the external magnetic field. In two particular cases, namely for *α* = 0 (i.e., in the CM geometry), and for *α* = π/2 (in the Faraday geometry) one obtains:

The total (elliptic) birefringence of the medium induced by an arbitrarily oriented magnetic field can be defined as [52]:

*n*and Δ

_{l}*n*are the linear and circular birefringences, i.e., the differences in the refractive index of the two wave components induced by the CM and the Faraday effect, respectively.

_{c}It is worth to recall that, in the Faraday geometry (i.e., for a beam propagating along the direction of the external magnetic field), the induced circular birefringence is proportional to the absolute value of the magnetic field (*b* ~|** H**|), while the linear birefringence caused by the CM effect (emerging when the direction of the light propagation is perpendicular to the magnetic field) is proportional to the square of the field,

*c*~|

**|**

*H*^{2}[52,53]. This implies a parabolic dependence between the birefringence coefficients for a given value of the magnetic field below saturation,

*c*~

*b*

^{2}, which offers the possibility of controlling the combined birefringences strength. In the general case of an oblique orientation of the magnetic field with respect to the propagation direction (assuming that the orientation is confined to a plane), i.e., with components of the magnetic field

*H*_{||}= |

**| cos**

*H**α*and

*H*_{⊥}= |

**| sin**

*H**α*(which are parallel and perpendicular to the direction of the light propagation, respectively), the total birefringence can be written as:

Following the considerations above, the strongest mixing between the polarization components is expected when the beam propagation direction is at an angle of 45° with respect to ** H**, which may significantly affect the collapse scenario for the bimodal beam. Moreover, both birefringence coefficients

*b*and

*c*, introduced in the theoretical model, are different from zero for an oblique direction of the magnetic field, whereas in the CM geometry

*b*= 0 and in the Faraday geometry

*c*= 0. However, in the vicinity of the collapse point, the paraxial approximation, which was assumed in the derivation of Eqs. (1), ceases to hold and the transversal components of the wave vector may appear even if they were previously absent. This means that in the CM geometry, a non-negligible contribution of the Faraday term (or equivalently, the CM term in the Faraday geometry) should be considered. As a result, more complex collapse dynamics may be expected in real conditions, as shown by our recent experimental data [61].

The magnetically-controlled linear and circular birefringences may be relatively easily induced in MO materials by a proper application of a magnetic field. For this purpose, magnetic crystals based on iron and rare-earth garnets may be used. In particular, iron garnets are cubic ferrimagnetic insulators described by the general formula, T_{3}Fe_{5}O_{12}, where T is a trivalent metal. The best-known representative of this class of materials is yttrium iron garnet (Y_{3}Fe_{5}O_{12}, abbreviated as YIG), which is highly transparent to near-infrared radiation, and is commonly used for the fabrication of optical isolators [62,63]. The nonlinear coefficients of YIG are expected to be relatively large (higher than in typical glasses), as suggested by the semi-empiric Miller’s rule, which states that a large refractive index implies the existence of a large nonlinear coefficient [3]. Undoped YIG at the wavelength of 1.2μm has typical values of *n*_{0} = 2.22, ∆*n _{F}* ≈1.6·10

^{−4}and ∆

*n*≈1·10

_{CM}^{−4}, at the saturation field of ~0.2T [53–55]. In such a case, the maximum values of the birefringence coefficients can be estimated as

*b*= 0.125,

*c*= 0.075 for a beam waist

*w*

_{0}= 5μm, in the Faraday and CM geometries, respectively. These values (and the ratio between them) can be tuned by adjusting the value and direction of the magnetic field, as well as by varying the spot-size and the wavelength of the injected beam (keeping also in mind that the circular magnetic birefringence strongly depends on the wavelength, while the linear magnetic birefringence is wavelength-independent for near-infrared radiation). Moreover, the range of the achievable birefringence parameters can be expanded using different magnetic materials. While in YIG only magnetic ions are ferric, magnetic rare-earth ions can be used in place of yttrium in order to improve the magnetic properties and thus increase the magneto-optic coefficients [53,63,64]. For example, the substitution of yttrium by cerium, resulting in Ce:YIG (with proper concentrations of Ce), may lead to values of the linear magnetic birefringence at least four times higher, and approximately six times larger for the circular birefringence, when compared to values observed in YIG at saturation [63,64].

## 3. Numerical results

#### 3.1. The input conditions

Generally, evolution equations in the form of Eqs. (1) support solutions that collapse after a finite propagation distance. In the case of MO bulk media, our objective is to understand whether the combination of the birefringence terms included in the considered model can essentially delay (or suppress) the collapse. Since the external magnetic field can change the total material birefringence, we expect it to affect the evolution of high-power beams in MO crystals with appropriate coefficients. In order to study the stabilizing effects of the birefringences, a detailed analysis of the dynamics of the collapse has been performed via the numerical solution of Eqs. (1), and, in parallel, of the variational equations, Eqs. (7). The propagation of beams with initially axisymmetric Gaussian profiles, described by Eq. (4), has been simulated. This choice of input was suggested by the fact that collapsing Gaussian beams are known to form the *Townes profile* that remains axisymmetric, irrespective of the initial spatial profile of the beam [8–10]. Moreover, the assumed axisymmetric structure of the vorticity-free solutions is not subject to azimuthal modulational instability [13,25]. The parameters for the input beam, Eq. (4), were taken as *ψ*_{0} = 0, *θ*_{0} = π/4, and *ξ*_{0} = 0. Thus, the profiles of the components at the input are given by:

*w*

_{0}and amplitude

*A*

_{0}.

The respective input total power (in physical units) is:

In order to avoid an artificial singularity at $r\to 0$ in Eqs. (1), the boundary conditions are taken as ${d{u}_{i}(r)/dr\text{\hspace{0.17em}}|}_{\text{\hspace{0.17em}}r=0}=0$ and ${\mathrm{lim}}_{\text{\hspace{0.17em}}r\to \infty}{u}_{i}(r)=0$, with the latter one implying that the beam is self-trapped (localized in the transverse plane). Moreover, the integration steps in both the transverse and the propagation directions, used in our numerical procedures, have been adjusted to ensure the conservation of the total power and of the Hamiltonian. As such, the very small deviations observed are contained within acceptable limits.

#### 3.2. Direct simulations (DS)

Firstly, numerical simulations of the evolution of the optical fields were performed via a direct solution of Eqs. (1), using an implicit axisymmetric central-difference Crank-Nicholson scheme [65]. We have checked that in most cases a resolution of $\Delta \tilde{r}$ = 10^{−3} over the transverse interval of $\tilde{r}$ = [0, 50] (and $\tilde{r}$ = [0, 100] for diffracting beams) together with a fixed propagation step-size of $\Delta \tilde{z}$ = 10^{−3}, allowed us to keep the relative deviation of the total power below 10^{−3} (with respect to the input value), over the considered distance (*z _{max}* = 100

*L*). We have also assumed that, for such a high-density meshing, a change in the critical power for the collapse caused by the numerical discretization [31] is negligible. In some particular cases (e.g., those in Fig. 2(c) and Fig. 3 ) the simulations were run over longer propagation distances, and the results were checked by re-running the simulations with a reduced mesh-grid-size, to ensure that the obtained results were independent from the chosen propagation and transverse step-sizes. Specifically, for the results represented by the peaks in Fig. 3(a1), 3(a2), a double-density mesh was used to make sure that none of observed effect was due to a coarse numerical scheme. Figure 1 presents the results obtained for two different input powers (here, the beam amplitude,

_{D}*A*

_{0}, is varied, while the width

*w*

_{0}is kept constant), as well as for different combinations of the magnetically-induced birefringences. In addition to the beam-intensity distributions (which are displayed in Figs. 1(a)–1(e) and are proportional to$\sum}_{i=R,L}{\left|\text{\hspace{0.05em}}{u}_{i}\text{\hspace{0.05em}}\right|}^{2$), the light propagation dynamics can also be represented by the normalized beam waist radius, $\tilde{R}(z)\equiv R(z)\text{\hspace{0.05em}\hspace{0.17em}}/\text{\hspace{0.17em}}R(z=0)$ [see Fig. 1(f)], where

*R*is defined as:

When the effects of the external magnetic field are neglected (i.e., *b* = *c* = 0), the beam polarization components are coupled only via the XPM term (which becomes significant only close to the collapse point). In such a case, the diffraction of a low-input power beam is observed as shown in Fig. 1(a). The spatial divergence of the beam may be reduced by increasing the optical power. For a given initial width, if the input amplitude exceeds a critical value (*A _{cr}* ≈1.124 for

*b*=

*c*= 0), the collapse occurs after a finite propagation distance [in particular, at

*z*≈12

*L*in Fig. 1(b)]. The numerically found critical amplitude,

_{D}*A*= 1.124, if rescaled back to physical units, corresponds to an optical power of

_{cr}*P*≈9.8MW (taking

*n*

_{0}= 2.22, λ

_{0}= 1.2μm and

*n*

_{2}= 10

^{−20}m

^{2}/W). As shown in Fig. 1(d), the introduction of the circular birefringence alone (i.e., b ≠ 0 while c = 0) does not change the collapse distance. However, the collapse can be accelerated (i.e., it can occur after a shorter propagation distance) by the introduction of the linear birefringence, c ≠ 0, as shown in Fig. 1(c). Contrarily, when both birefringences are present, the collapse distance may be extended, as shown in Fig. 1(e). The latter outcome takes place only for a particular combination of the birefringence coefficients [see Fig. 2(c)], resulting from a specific interplay between the amplitude and phase mixing of the two polarization components. Thus, taking into account that the coefficients introduced in the model depend on the value of the external magnetic field, as described in section 2.3, the distance for the collapse may be magnetically tuned, allowing the possibility of controlling the self-focusing beams in bulk MO media.

The effect of the combined birefringences for a high-intensity beam with a fixed amplitude (*A*_{0} = 1.135, which is above the critical value) is schematically presented in Fig. 2, where the total propagation length is limited to *z _{max}* = 100

*L*. The figure demonstrates the tunability of the collapse in the (

_{D}*b*,

*c*)-plane. In particular, Fig. 2(a) shows the normalized final width of the beam (obtained by fitting the spatial profile of one of the beam components to a Gaussian function) at

*z*=

*z*. If the collapse occurs at

_{max}*z*≤

*z*, the final width of the beam is set to be

_{max}*w*= 0 [the grey area in Fig. 2(a)]. The corresponding normalized collapse distance,

*z*/

*L*, is shown in Fig. 2(b). The results are displayed for the

_{D}*u*component, while the dynamics of

_{R}*u*is essentially the same in the entire range of the considered parameters.

_{L}The black regions in Fig. 2(a) [and also in Fig. 2(d)] indicate the output widths that are close to the initial value, i.e., *w*/*w*_{0} ≈1, thus defining the range of parameters where the propagating beams are potentially stable. As illustrated in Fig. 2(b), at a given value of the input power, the amplitude mixing between the polarization components dominates for *c* > *b*, which usually results in the acceleration of the collapse, while a prominent phase mixing (for *b* > *c*) may lead to the suppression of the collapse. Furthermore, distinctive regions in the chart can be singled out, as labeled in Fig. 2(a). Of particular importance are three characteristic (critical) values of the circular-birefringence coefficient *b* (for a fixed value of *c*): these are marked as *b _{A}*,

*b*and

_{B}*b*in the figures. Two of them (

_{C}*b*and

_{A}*b*) correspond to the transition between the collapsing and diffracting solutions. The third critical point,

_{C}*b*, indicates parameters for which the beam divergence may be reduced, leading to a longer stable propagation length

_{B}*L*. Eventually, as one can see in Fig. 2(c), the length for which a nearly stable propagation takes place (before the beam eventually collapses or diffracts) may be significantly extended by precisely choosing the proper combination of the birefringence coefficients (

*b*,

*c*), and in particular, it can be increased to up to about 300 diffraction lengths in the specific case presented here.

The dynamics in the proximity to the critical points is illustrated in Fig. 3. It is apparent that the mechanism of the collapse delay at the diffraction-collapse border (i.e., for the cases of *b _{A}* and

*b*) is independent of the linear birefringence coefficient

_{C}*c*. The obtained numerical results at these boundaries were determined to be well represented by power-functions to approximate the collapse length,

*L*:

In particular, for the cases of *c* = 0.055 and *c* = 0.07 [namely those for *A* = 1.135 presented in Fig. (3)], the critical exponents *α _{A,C}* were found to be around 0.05, indicating rapid changes (increasing with input optical power) of the collapse length

*L*close to the collapse-diffraction transition point. Thus, a strong delay of the collapse is possible for the parameters that are close enough to the critical values determined by

*b*and

_{A}*b*. For a given linear birefringence coefficient

_{C}*c*, a prolongation of the stable propagation of the beam is also possible if the value of the co-acting circular birefringence coefficient

*b*is approaching the above-mentioned third critical value,

*b*. In such cases, the range of suitable values for the propagation stabilization,

_{B}*∆b*, broadens as a function of the linear birefringence coefficient

_{B}*c*[see Figs. 3(a1) and 3(a2)]. The normalized beam waist radius as a function of the propagation distance for the birefringence coefficients approaching these critical values was also numerically determined, and is presented in Fig. 3(b1), 3(b2).

Moreover, the critical circular birefringence coefficients *b _{B}* feature a nearly linear dependence on

*c*, as shown in Fig. 4(a) . This plot allows one to determinate the required combination of the MO parameters for which the beam divergence may be reduced [examples are shown in Fig. 4(b1), 4(b2)], and, specifically, to identify situations [black squares in Fig. 4(a)] when the nearly stable propagation (before the eventual diffraction) can be extended. The quasi-stabilization line obtained for

*A*

_{0}= 1.135 is shown in Fig. 4(a), while the evident prolongation of the stable propagation [similar to that shown in Figs. 3(b1), 3(b2)] is observed only for a limited range of the combined birefringences, marked with the black squares in Fig. 4(a).

#### 3.3. The variational approximation

Variational equations [Eqs. (7)] have been solved numerically using an explicit, variable-step-size Runge-Kutta (RK) scheme, implemented with the use of the Dormand-Prince pair method. In particular, the truncation error was estimated by comparing the results produced using fourth- and fifth-order RK formulas, in order to dynamically adjust the step-size [65,66]. The input Gaussian profile defined by Eq. (12) was used as the initial condition. We have concluded that the maximum step-size allowed by the numerical routine, $\Delta \tilde{z}$ = 0.01, ensures the conservation of the Hamiltonian [given by Eq. (8)] with relative differences below 10^{−11}. In the case of collapse, the beam width varies so fast that even a drastic reduction of the step-size does not secure the conservation of the Hamiltonian. In particular, we have defined the collapse point as the propagation distance for which the numerical procedure fails to stay within the predetermined margins, fixed to 10^{−3} and 10^{−6} for the relative and absolute error tolerances, respectively, without reducing the step-size below the smallest value allowed by the numerical routine, $\Delta \tilde{z}$ ~10^{−14}.

The simulations we have performed demonstrate that in the absence of a magnetic field, the collapse takes place when the beam power reaches the critical value, ${\left({\tilde{N}}_{cr}\right)}_{VA}=2/3$ (this well-known [59] critical value is immediately predicted by Eq. (7a) with ${\mathrm{sin}}^{2}\left(2\theta \right)=1$, which implies a linear polarization, as considered above). When comparing to the results obtained from the direct simulations of Eqs. (1), with *b* = *c* = 0 and $\theta =\text{\pi}/4$, we find that the critical power is larger using the VA, with their ratio being equal to${({N}_{cr})}_{DS}/{({N}_{cr})}_{VA}\approx 0.947$. This implies that the input beam amplitude, *A*_{0} = 1.135, used for the PDE simulations, corresponds to the initial condition of ${\tilde{N}}_{0}=0.68$ for the VA.

The numerical results obtained for this value of ${\tilde{N}}_{0}$ are summarized in Fig. 2(d) and Fig. 4(a). It is evident that both the PDE- and VA-based simulations yield quite similar results, and nearly identical qualitative predictions regarding the possibility to control the beam self-focusing and collapse by means of the magnetically-induced birefringences. In particular, the range of combined birefringence coefficients for which the collapse may be prevented (for the analyzed value of power) is almost the same, as can be seen from Figs. 2(a) and 2(d). Moreover, the slopes of the (quasi-) stabilization lines, shown in Fig. 4(a), are nearly identical. Hence, the VA allows one to capture the main physical effects and predict the behavior of the beam propagating in the MO medium, with much shorter computational times required for the simulations.

Nevertheless, certain differences between the results obtained by the PDE and VA approaches are observed for parameters close to their critical values, and they get more pronounced when the values of the birefringence coefficients increase. In particular, this pertains to the width of the stabilization region, ∆*b _{B}*, which is easily observed in the case of the VA; see the black region in Fig. 2(d). Moreover, as it is presented in Fig. 5(a)
, several sharp maxima of the propagation length within the stabilization region are obtained for critical values

*b*, in contrast to the DS results [Fig. 3(a1), 3(a2)]. The width of the stabilization range, ∆

_{Bi}*b*, as well as the number of peak values corresponding to the stable-propagation lengths, increases with the value of the linear birefringence coefficient

_{B}*c*(not shown here). Moreover, while the solutions found for values of

*b*taken in the vicinity of

*b*exhibit oscillations along the propagation distance with a constant period (as it was the case for the DS), irregular oscillations of the beam width for

_{A}*b*close to

*b*and

_{Bi}*b*are observed [see Fig. 5(b)]. Specifically, solutions around

_{C}*b*are prone to instabilities, and the transition between the diffraction and collapse is not as sharp as it was for the DS. In such a case collapse is observed for

_{C}*b*≥

*b*but the width of the beam exceed significantly its initial value before the collapse point [as seen on the bottom graph in Fig. 5(b) and as indicated by the solid dots in Fig. 5(a)]. However, it is still possible to estimate the critical values of the circular birefringence coefficient, and to fit it to the power-function approximation given by Eq. (20) (for instance,

_{C}*α*≈0.07 and

_{A}*α*≈0.12 in the present case).

_{C}One of the reasons for the different results produced by the VA and the PDE simulations (e.g., the value retrieved for the collapse distance) is due to the fact that the VA assumes that the spatial-beam profile remains Gaussian, while the direct simulations allow for a change of the beam shape along propagation (in the full simulations, the shape indeed deviates from a perfect Gaussian).

## 4. Conclusions

In conclusion, we have reported the results of the numerical analyses performed for the propagation of nonlinear bimodal beams in bulk Kerr-type magneto-optical media. In particular, the possibility of modifying the collapse dynamics by means of an external magnetic field was demonstrated. Direct PDE simulations, as well as investigations based on the variational approximation have been carried out, yielding consistent results. We have also demonstrated that the acceleration, the delay, and even the arrest of the collapse may be achieved via a proper combination of the linear and the circular birefringences, induced by the magnetic field through the Cotton-Mouton and Faraday effects, respectively. The ability of controlling the catastrophic self-focusing of optical beams is of great importance for practical applications and can be readily accomplished in available magnetic materials, such as rare-earth garnets (e.g., see [61]). In addition, our analysis suggests new possibilities for the implementation of MO-driven effects in nonlinear optics.

## Acknowledgments

Katarzyna Rutkowska would like to acknowledge a Marie Curie fellowship (MOIF-CT-2006-039600). Roberto Morandotti would like to thank the NSERC (Natural Science and Engineering Research Council) Strategic Grant Program, as well as TeraXion Ltd and O E Land Ltd for support of this work.

## References and links

**1. **C. Sulem, and P. L. Sulem, *The nonlinear Schrödinger equation: self-focusing and wave collapse* (Springer, New York 1999).

**2. **E. Infeld, and G. Rowlands, *Nonlinear Waves, Solitons and Chaos* (Cambridge Univ. Press, Cambridge 2000).

**3. **R. W. Boyd, *Nonlinear Optics* (Academic Press, San Diego 2008).

**4. **B. A. Malomed, “Nonlinear Schrödinger equations,” in *Encyclopedia of Nonlinear Science*, A. Scott, ed. (Routledge, New York 2005).

**5. **Y. S. Kivshar, and P. G. Agrawal, *Optical Solitons: From Fibers to Photonic Crystal* (Academic Press, San Diego 2003).

**6. **S. Trillo, and W. Torruellas, *Spatial Solitons* (Springer-Verlag, Berlin 2001).

**7. **G. I. Stegeman and M. Segev, “Optical Spatial Solitons and Their Interactions: Universality and Diversity,” Science **286**(5444), 1518–1523 (1999). [CrossRef] [PubMed]

**8. **R. Y. Chiao, E. Garmire, and C. H. Townes, “Self-Trapping of Optical Beams,” Phys. Rev. Lett. **13**(15), 479–482 (1964). [CrossRef]

**9. **K. D. Moll, A. L. Gaeta, and G. Fibich, “Self-similar optical wave collapse: observation of the Townes profile,” Phys. Rev. Lett. **90**(20), 203902 (2003). [CrossRef] [PubMed]

**10. **G. Fibich and A. L. Gaeta, “Critical power for self-focusing in bulk media and in hollow waveguides,” Opt. Lett. **25**(5), 335–337 (2000). [CrossRef] [PubMed]

**11. **M. Weinstein, “Nonlinear Schrödinger Equations and Sharp Interpolation Estimates,” Commun. Math. Phys. **87**(4), 567–576 (1983). [CrossRef]

**12. **T. Lahaye, J. Metz, B. Fröhlich, T. Koch, M. Meister, A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi, and M. Ueda, “d-wave collapse and explosion of a dipolar bose-einstein condensate,” Phys. Rev. Lett. **101**(8), 080401 (2008). [CrossRef] [PubMed]

**13. **T. Koch, T. Lahaye, J. Metz, B. Fröhlich, A. Griesmaier, and T. Pfau, “Stabilization of a purely dipolar quantum gas against collapse,” Nat. Phys. **4**(3), 218–222 (2008). [CrossRef]

**14. **P. L. Kelley, “Self-focusing of optical beams,” Phys. Rev. Lett. **15**(26), 1005–1008 (1965). [CrossRef]

**15. **L. Bergé, “Wave collapse in physics: principles and applications to light and plasma waves,” Phys. Rep. **303**(5-6), 259–370 (1998). [CrossRef]

**16. **T. D. Grow, A. A. Ishaaya, L. T. Vuong, and A. L. Gaeta, “Collapse and stability of necklace beams in Kerr media,” Phys. Rev. Lett. **99**(13), 133902 (2007). [CrossRef] [PubMed]

**17. **D. P. Lathrop, B. W. Zeff, B. Kleber, and J. Fineberg, “Singularity dynamics in curvature collapse and jet eruption on a fluid surface,” Nature **403**(6768), 401–404 (2000). [CrossRef] [PubMed]

**18. **P. A. Robinson, “Nonlinear wave collapse and strong turbulence,” Rev. Mod. Phys. **69**(2), 507–574 (1997). [CrossRef]

**19. **Y. R. Shen, “Recent advances in nonlinear optics,” Rev. Mod. Phys. **48**(1), 1–32 (1976). [CrossRef]

**20. **E. A. Kuznetsov, J. J. Rasmussen, K. Rypdal, and S. K. Turitsyn, “Shaper criteria for the wave collapse,” Physica D **87**(1-4), 273–284 (1995). [CrossRef]

**21. **Y. Silberberg, “Collapse of optical pulses,” Opt. Lett. **15**(22), 1282–1284 (1990). [CrossRef] [PubMed]

**22. **F. Wise and P. di Trapani, “The Hunt for Light Bullets – Spatiotemporal Solitons,” Opt. Photonics News **13**(2), 29 (2002).

**23. **B. A. Malomed, D. Mihalache, F. Wise, and L. Torner, “Spatiotemporal optical solitons,” J. Opt. B: Quant. Semicl. Opt. **7**(5), R53–R72 (2005). [CrossRef]

**24. **A. Chong, W. H. Renninger, D. N. Christodoulides, and F. W. Wise, “Airy-Bessel wave packets as versatile linear light bullets,” Nat. Photonics **4**(2), 103–106 (2010). [CrossRef]

**25. **Y. S. Kivshar and D. E. Pelinovsky, “Self-focusing and transverse instabilities of solitary waves,” Phys. Rep. **331**(4), 117–195 (2000). [CrossRef]

**26. **A. L. Gaeta, “Catastrophic collapse of ultrashort pulses,” Phys. Rev. Lett. **84**(16), 3582–3585 (2000). [CrossRef] [PubMed]

**27. **G. Fibich and G. C. Papanicolaou, “Self-focusing in the presence of small time dispersion and nonparaxiality,” Opt. Lett. **22**(18), 1379–1381 (1997). [CrossRef] [PubMed]

**28. **D. Cheskis, S. Bar-Ad, R. Morandotti, J. S. Aitchison, H. S. Eisenberg, Y. Silberberg, and D. Ross, “Strong spatiotemporal localization in a silica nonlinear waveguide array,” Phys. Rev. Lett. **91**(22), 223901 (2003). [CrossRef] [PubMed]

**29. **N. Akhmediev, A. Ankiewicz, and J. M. Soto-Crespo, “Does the nonlinear Schrödinger equation correctly describe beam propagation?” Opt. Lett. **18**(6), 411–413 (1993). [CrossRef] [PubMed]

**30. **G. Fibich and B. Ilan, “Optical light bullets in a pure Kerr medium,” Opt. Lett. **29**(8), 887–889 (2004). [CrossRef] [PubMed]

**31. **O. Bang, J. J. Rasmussen, and P. L. Christiansen, “Subcritical localization in the discrete nonlinear Schrödinger equation with arbitrary power nonlinearity,” Nonlinearity **7**(1), 205–218 (1994). [CrossRef]

**32. **A. Szameit, J. Burghoff, T. Pertsch, S. Nolte, A. Tünnermann, and F. Lederer, “Two-dimensional soliton in cubic fs laser written waveguide arrays in fused silica,” Opt. Express **14**(13), 6055–6062 (2006). [CrossRef] [PubMed]

**33. **O. Bang, D. Edmundson, and W. Królikowski, “Collapse of incoherent light beams in inertial bulk Kerr media,” Phys. Rev. Lett. **83**(26), 5479–5482 (1999). [CrossRef]

**34. **A. S. Desyatnikov, D. Buccoliero, M. R. Dennis, and Y. S. Kivshar, “Suppression of collapse for spiraling elliptic solitons,” Phys. Rev. Lett. **104**(5), 053902 (2010). [CrossRef] [PubMed]

**35. **L. Bergé, O. Bang, J. J. Rasmussen, and V. K. Mezentsev, “Self-focusing and solitonlike structures in materials with competing quadratic and cubic nonlinearities,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics **55**(3), 3555–3570 (1997). [CrossRef]

**36. **V. Skarka, V. I. Berezhiani, and R. Miklaszewski, “Spatiotemporal soliton propagation in saturating nonlinear optical media,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics **56**(1), 1080–1087 (1997). [CrossRef]

**37. **S. Gatz and J. Herrmann, “Propagation of optical beams and the properties of two-dimensional spatial solitons in media with a local saturable nonlinear refractive index,” J. Opt. Soc. Am. B **14**(7), 1795 (1997). [CrossRef]

**38. **O. Bang, W. Krolikowski, J. Wyller, and J. J. Rasmussen, “Collapse arrest and soliton stabilization in nonlocal nonlinear media,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. **66**(4), 046619 (2002). [CrossRef] [PubMed]

**39. **M. Peccianti, K. A. Brzdkiewicz, and G. Assanto, “Nonlocal spatial soliton interactions in nematic liquid crystals,” Opt. Lett. **27**(16), 1460–1462 (2002). [CrossRef] [PubMed]

**40. **P. Pedri and L. Santos, “Two-dimensional bright solitons in dipolar Bose-Einstein condensates,” Phys. Rev. Lett. **95**(20), 200404 (2005). [CrossRef] [PubMed]

**41. **I. Towers and B. A. Malomed, “Stable (2+1)-dimensional solitons in a layered medium with sign-alternating Kerr nonlinearity,” J. Opt. Soc. Am. B **19**(3), 537 (2002). [CrossRef]

**42. **F. Kh. Abdullaev, J. G. Caputo, R. A. Kraenkel, and B. A. Malomed, “Controlling collapse in Bose-Einstein condensation by temporal modulation of the scattering length,” Phys. Rev. A **67**(1), 013605 (2003). [CrossRef]

**43. **H. Saito and M. Ueda, “Dynamically stabilized bright solitons in a two-dimensional bose-einstein condensate,” Phys. Rev. Lett. **90**(4), 040403 (2003). [CrossRef] [PubMed]

**44. **S. Eisenmann, E. Louzon, Y. Katzir, T. Palchan, A. Zigler, Y. Sivan, and G. Fibich, “Control of the filamentation distance and pattern in long-range atmospheric propagation,” Opt. Express **15**(6), 2779–2784 (2007). [CrossRef] [PubMed]

**45. **J. Yang and Z. H. Musslimani, “Fundamental and vortex solitons in a two-dimensional optical lattice,” Opt. Lett. **28**(21), 2094–2096 (2003). [CrossRef] [PubMed]

**46. **C. R. Menyuk, “Stability of solitons in birefringent optical fibers. I: equal propagation amplitudes,” Opt. Lett. **12**(8), 614–616 (1987). [CrossRef] [PubMed]

**47. **S. Trillo, S. Wabnitz, E. M. Wright, and G. I. Stegeman, “Polarized soliton instability and branching in birefringent fibers,” Opt. Commun. **70**(2), 166–172 (1989). [CrossRef]

**48. **B. A. Malomed, “Polarization dynamics and interactions of solitons in a birefringent optical fiber,” Phys. Rev. A **43**(1), 410–423 (1991). [CrossRef] [PubMed]

**49. **Y. Barad and Y. Silberberg, “Polarization Evolution and Polarization Instability of Solitons in a Birefringent Optical Fiber,” Phys. Rev. Lett. **78**(17), 3290–3293 (1997). [CrossRef]

**50. **L. Bergé, O. Bang, and W. Królikowski, “Influence of four-wave mixing and walk-Off on the self-focusing of coupled waves,” Phys. Rev. Lett. **84**(15), 3302–3305 (2000). [CrossRef] [PubMed]

**51. **O. Bang, L. Bergé, and J. J. Rasmussen, “Fusion, collapse, and stationary bound states of incoherently coupled waves in bulk cubic media,” Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics **59**(4), 4600–4613 (1999). [CrossRef]

**52. **R. Kurzynowski and W. A. Woźniak, “Superposition rule for the magneto-optic effects in isotropic media,” Optik (Stuttg.) **115**(10), 473–475 (2004). [CrossRef]

**53. **A. K. Zvezdin, and V. A. Kotov, *Modern Magnetooptics and Magnetooptical Materials* (Taylor & Francis Group, New York 1997).

**54. **J. F. Dillon Jr, J. P. Remeika, and C. R. Staton, “Linear Magnetic Birefringence in the Ferrimagnetic Garnets,” J. Appl. Phys. **41**(11), 4613 (1970). [CrossRef]

**55. **G. B. Scott, D. E. Lacklison, H. I. Ralph, and J. L. Page, “Magnetic circular dichroism and Faraday rotation spectra of Y_{3}Fe_{5}O_{12},” Phys. Rev. B **12**(7), 2562 (1975). [CrossRef]

**56. **G. A. Smolenskii, R. V. Pisarev, and I. G. Sinii, “Birefringence of light in magnetically ordered crystals,” Usp. Fiziol. Nauk **116**, 231 (1975). [CrossRef]

**57. **J. M. Dudley, G. Genty, and S. Coen, “Supercontinuum generation in photonic crystal fiber,” Rev. Mod. Phys. **78**(4), 1135–1184 (2006). [CrossRef]

**58. **R. J. Ballagh, K. Burnett, and T. F. Scott, “Theory of an Output Coupler for Bose-Einstein Condensated Atoms,” Phys. Rev. Lett. **78**(9), 1607–1611 (1997). [CrossRef]

**59. **B. A. Malomed, “Variational methods in nonlinear fiber optics and related fields,” Progr. Opt. **43**, 71 (2002). [CrossRef]

**60. **M. Desaix, D. Anderson, and M. Lisak, “Variational approach to collapse of the optical pulses,” J. Opt. Soc. Am. B **8**(10), 2082 (1991). [CrossRef]

**61. **Y. Linzon, K. A. Rutkowska, B. A. Malomed, and R. Morandotti, “Magneto-optical control of light collapse in bulk Kerr media,” Phys. Rev. Lett. **103**(5), 053902 (2009). [CrossRef] [PubMed]

**62. **J.-M. Liu, *Photonic Devices*, Cambridge Univ. Press, New York 2005.

**63. **O. Kamada, T. Nakaya, and S. Higuchi, “Magnetic field optical sensors using Ce:YIG single crystals as a Faraday element,” Sens. Actuators A Phys. **119**(2), 345–348 (2005). [CrossRef]

**64. **M. C. Sekhar, J.-Y. Hwang, M. Ferrera, Y. Linzon, L. Razzari, C. Harnagea, M. Zaezjev, A. Pignolet, and R. Morandotti, “Strong enhancement of the Faraday rotation in Ce: and Bi: co-modified epitaxial iron garnet thin films,” Appl. Phys. Lett. **94**(18), 181916 (2009). [CrossRef]

**65. **U. M. Ascher, *Numerical methods for evolutionary differential equations* (SIAM, Philadelphia 2008).

**66. **J. C. Butcher, *Numerical Methods for Ordinary Differential Equations* (Wiley & Sons, West Sussex 2008).