1. Introduction

In this paper, the rotationally symmetric finite-difference time-domain (FDTD) formulation of Maloney, Smith and Scott’s method [1], was integrated with Greene and Taflove’s FDTD general vector auxiliary differential equation (GVADE) method [2–6], allowing for a radially polarized 3D cylindrical rotationally symmetric (CRS) electromagnetic wave to be simulated in nonlinear dispersive material using a 2D numerical grid. An example of a 3D CRS soliton is shown in Fig. 1. The FDTD GVADE method for modeling electromagnetic vector wave propagation in 3D nonlinear material requires a 3D grid to simulate a 3D electromagnetics wave, assuming the linear behavior is accurately represented by Lorentz polarization and the nonlinearity is accurately represented by Kerr and Raman expressions incorporating third order susceptibility. In comparison, [1] uses FDTD to simulate a 3D CRS source and a 3D CRS propagating electromagnetic wave in a linear material with a 2D grid by taking advantage of rotational symmetry. The integration of these methods allows the simulation of 3D CRS electromagnetic waves in nonlinear material using 2D computational resources, which is a significant computational improvement in speed and memory. To better understand the 3D CRS behavior, below we present the 2D theory and results alongside the 3D CRS theory and results.

 

Fig. 1. An example of a cylindrical rotationally symmetric spatial soliton. This figure is a plot of the magnitude of the electric field $|{\textbf E} |$, showing the 3D aspect of the wave as it propagates in the positive z direction. The simulation results used for this figure are the same as Fig. 4(c), where the fields are only plotted in the y = 0 plane for comparison with 2D simulation results. Note: The slice of the fields removed from $\phi \mathrm{\ =\ \ 0^\circ \ \; to\;\ }\phi \mathrm{\ =\ \ 90^\circ } $ are not plotted for visual purposes, but the actual wave would be a complete cylinder without the slice removed.

Download Full Size | PDF

2. Theory

The method for solving 3D transverse magnetic to z ($\textrm{T}{\textrm{M}_\textrm{z}}$) CRS equations for the electromagnetic fields propagating in a nonlinear medium at optical frequencies, involves combining the 3D CRS FDTD update equations and rotational symmetry assumptions of [1], which builds on [7,8], with the 2D $\textrm{T}{\textrm{M}_\textrm{z}}$ version of the FDTD GVADE method [2,3,5], which generalized past FDTD auxiliary differential equation (ADE) methods [9–11]. First the FDTD GVADE method numerically solves Maxwell’s curl equations, Eq. (1) and Eq. (2). In (1) and (2), $\textbf{H}$ and E are the magnetic and electric field vectors and ${\textbf{J}_{\textrm{Kerr}}}$, ${\textbf{J}_{\textrm{Raman}}}$ and ${\textbf{J}_{\textrm{Lorentz,l}}}$ are polarization current vectors incorporated in Maxwell’s equations according to,

Starting with the magnetic field at time step n + 1/2 from Ampere’s law Eq. (1), the electric fields at time-step n + 1 is found through an application of the multi-dimensional Newton-Raphson method that includes numerical expressions for the material parameters, linear and nonlinear, accounted for by the polarization currents. Next, with Faraday’s law, shown in Eq. (2), expressed as a finite difference update equation, the new electric fields are used to calculate the magnetic field advanced one-half time-step forward in time. When solving Ampere’s Law for the electric field using the multidimensional Newton-Raphson method, the FDTD GVADE method first moves the curl term in Ampere’s Law to the right-hand side of the equation,

Since the 2D form is the most relevant for this paper, a 2D $\textrm{T}{\textrm{M}_\textrm{z}}$ version of this equation is shown in Eq. (5). When the field components do not depend on the y coordinate, this causes the $\frac{\partial }{{\partial \textrm{y}}}$ terms from Maxwell’s equations to be zero, leaving only the ${\textrm{H}_\textrm{y}}$, ${\textrm{E}_x}$ and ${E_z}$ fields,

This equation can be re-written in a form discretized with time-steps (n) using finite differences, which can be solved by the Newton-Raphson method, whose goal is to minimize X, Y and Z, ideally making them zero, to thereby determine the electric field ${\textbf{E}^{{n + 1}}}$,

After the electric fields ${\textbf{E}^{{n + 1}}}$ are determined using the Newton-Raphson method, the FDTD GVADE method proceeds to the next timestep (n = n + 1), where the ${\textbf{E}^{{n + 1}}}$ fields are now termed the new ${\textbf{E}^{n}}$ fields. Next, the ${\textbf{E}^{n}}$ fields are plugged into the standard FDTD magnetic field update equations based on Faradays Law shown in Eq. (2) to determine the new ${\textbf{H}^{{n + }\frac{\textrm{1}}{\textrm{2}}}}$, which is used to repeat the process to find the new ${\textbf{E}^{{n + 1}}}$ [2–6].

For the case of the 3D $\textrm{T}{\textrm{M}_\textrm{z}}$ CRS waves, the only changes required are converting the coordinate system from Cartesian to cylindrical, which changes the $\nabla \mathrm{\ \times }{\textbf{H}^{{n + }\frac{\textrm{1}}{\textrm{2}}}}$ term, thereby enforcing the cylindrically symmetric boundary conditions at r = 0, and implementing a 3D rotationally symmetric source/excitation. Re-writing Eq. (5) in cylindrical coordinates, when the field components do not depend on the $\phi $ coordinate, this causes the $\frac{\partial }{{\partial \phi }}$ terms from Maxwell’s equations to be zero, leaving only ${\textrm{H}_\phi }$, ${\textrm{E}_r}$ and ${E_z}$ fields,

In the equation above R represents the radial component of the electric field vector. The following sub-sections further explain the theory of the 3D $\textrm{T}{\textrm{M}_\textrm{z}}$ CRS FDTD GVADE method, and compare it with the 2D Cartesian version.

2.1 Cartesian 2D waves

For comparison’s sake with 3D CRS, the 2D $\textrm{T}{\textrm{M}_\textrm{z}}$ Cartesian equations are shown below:

The simplified Ampere’s law Eqs. (10) and (11) are solved using the 2D Newton-Raphson method for ${{E}_x}$ and ${{E}_{z}}$ after being re-written as FDTD finite difference update equations,

The remaining Eq. (9) can be solved for the term including ${{H}_y}$, and then resolved for ${H}_y^{\; n + 1/2}$, after re-writing it as a standard FDTD update equation.

2.2 Cylindrical rotationally symmetric 3D waves

The similar 3D $\textrm{T}{\textrm{M}_\textrm{z}}$ CRS equations are shown below:

The simplified Ampere’s law Eqs. (15) and (16) are solved using the 2D Newton-Raphson method for ${{E}_{r}}$ and ${{E}_{z}}$ after being re-written as FDTD finite difference update equations,

Using the product rule, the last equation can be re-written in another form:

The remaining Eq. (14) can be solved for the term including ${H_\phi }$, and then resolved for $H_\phi ^{\; n + 1/2}$, after re-writing it as a standard FDTD update equation.

2.3 Comparing 2D Cartesian and 3D rotationally symmetric waves

The 2D Cartesian and 3D CRS equations are quite similar in form, and it can be observed that they are equivalent in the y = 0 plane for ${x\; } \ge {\; 0}$ ($\phi \mathrm{\;\ =\ \;\ 0^\circ }$), resulting in ${r\; = \; }\sqrt {{{x}^{2}}{ + }{{y}^{2}}} {\; = \; }\sqrt {{{x}^{2}}} { = \; |x|}$ and z = z, with $\phi $ and y both pointing “out of the page”. The only difference in form is the “extra ${{H}_\phi }$/r term” in the ${({\nabla \times {\textbf{H}}} )_{{z,3D}}}$ equation compared with the ${({\nabla \times {\textbf{H}}} )_{{z,2D}}}$ equation as shown below:

Re-writing the equations in the y = 0 plane,

In the y = 0 plane, if symmetry about the line x = 0, the z-axis, is taken into account, even symmetry for ${{E}_{z}}$, odd symmetry for ${{H}_{y}}$ (${{H}_\phi }$) and ${{E}_{x}}$ (${{E}_{r}}$), this allows us to cut the computational space in half in the x direction, resulting in the advantage that only the ${x\; } \ge {\; 0}$ region ($\phi \mathrm{\;\ =\ \;\ 0^\circ }$) needs to be analyzed, assuming the source is centered around r = 0. The even and odd symmetry about the z-axis is the result of rotational symmetry for the 3D CRS case. Additionally, for both the 2D Cartesian and 3D CRS cases, the odd function nature about the z-axis of the FDTD GVADE Ampere’s law terms, $\nabla {\ \times \textbf{H}}$, $\frac{{\partial {\textbf{E}}}}{{\partial {t}}}$, ${{\textbf{J}}_{{\textrm{Kerr}}}}$, ${{\textbf{J}}_{{\textrm{Raman}}}}$ and ${{\textbf{J}}_{{\textrm{Lorentz},l}}}$ from Eq. (6) and Eq. (8), results in odd symmetry of the ${{E}_{x}}$ and ${{E}_{r}}$ waveforms since ${({\nabla \times {\textbf{H}}} )_{x{,2D}}}$ and ${({\nabla \times {\textbf{H}}} )_{r{,3D}}}$ are odd waveforms given an odd ${{H}_{y}}$ (${{{H}}_\phi }$) excitation/source waveform, while for an odd ${{{H}}_{y}}$ (${{H}_\phi }$) excitation/source waveform ${({\nabla \times {\textbf{H}}} )_{{z,2D}}}$ and ${({\nabla \times {\textbf{H}}} )_{{z,3D}}}$ are even waveforms resulting in even symmetry for the ${{E}_{z}}$ waveform. The symmetry is maintained as the waveforms propagate. This can be shown to be true using the known properties of adding, subtracting, and multiplying even and odd functions. These features can be seen by looking at Figs. 2, 8,  9, and 10. This results in the following equations in the y = 0 and ${x\; } \ge {\; 0}$ region ($\phi \mathrm{\;\ =\ \;\ 0^\circ }$), summarizing the validity of comparing the 2D and 3D CRS results in this region:

In the Results section, the 2D Cartesian simulation will be compared with the 3D CRS simulation results in the y = 0 plane where ${x\; } \ge {\; 0}$ ($\phi \mathrm{\;\ =\ \;\ 0^\circ }$), to show the effect of the added ${H_\phi }/r$ term in the 3D CRS simulations.

It is interesting to note that when r becomes large, the ${{H}_\phi }/{r}$ term becomes smaller, causing ${({\nabla \times {\textbf{H}}} )_{{\textrm{z},\; \textrm{3D}}}}$ to approach becoming equal to the ${({\nabla \times {\textbf{H}}} )_{{\textrm{z},\; \textrm{2D}}}}$ term in the y = 0 plane,

This implies that for ${{H}_\phi }$ distributions located at larger r values, the 3D behavior and the 2D behavior should be similar in the y = 0 plane. Some example cases of ${({\nabla \times {\textbf{H}}} )_{{\textrm{z},\; \textrm{3D}}}}$ becoming progressively closer to ${({\nabla \times {\textbf{H}}} )_{{\textrm{z},\; \textrm{2D},\; x\; } \ge {\; 0,\; {y}\; = \; 0}}}$ as ${{H}_\phi }/{r}$ decreases are shown for “source spacings” = 7.5, 10 and 15 in the Results section. The 3D behavior progressively approaches the known behavior for two antiphase, a phase difference of π, (1 + 1)D, termed “2D” for this paper, solitons with larger hyperbolic secant separation distances. The term “source spacing” is defined in the Excitation Waveform section.

Please note that the choice of the y = 0 plane, to compare the 2D and 3D CRS equations, was for convenience. Due to the rotational symmetry, the 3D CRS fields do not vary based on the choice of plane with cylindrical/azimuthal angle $\phi $ through the z-axis by which they are sampled, and Eqs. (14–19) do not lose their generality since the fields are merely being sampled or measured in this plane at this angle, not changed or simplified. Alternatively, as an example, the 2D Cartesian case could have been chosen to have no field variation in the x-direction instead of the y-direction, resulting in ${{E}_y}$, ${{E}_{z}}$ and ${{H}_{x}}$ fields with equations similar to Eqs. (9–11), allowing comparison with 3D CRS in the x = 0 plane instead of y = 0, where a similar set of equations to Eqs. (22–24) could be derived.

2.4 Cylindrical rotationally symmetric 3D wave boundary conditions at r = 0

By observation, rotational symmetry of ${{H}_\phi }$ and ${{E}_{r}}$ implies these field components must be continuous, transitioning from the $\phi \mathrm{\;\ =\ \;\ 0^\circ }$ plane to the $\phi \mathrm{\;\ =\ \;\ 180^\circ }$ plane at r = 0. As a result, ${{H}_\phi }{\; = \; 0}$ and ${{E}_r}{ = \; 0}$ at r = 0. In other words, ${{H}_\phi }$ and ${{E}_{r}}$ have “odd symmetry” across r = 0 when moving $\mathrm{180^\circ }$ from the starting angle $\phi $. For example, in the y = 0 plane, ${{H}_\phi }$ and ${{E}_{r}}$.have odd symmetry in x about x = 0 when a CRS source is centered at r = 0.

It can by shown by rotational symmetry and L’Hopital’s rule that at r = 0 the boundary conditions are:

2.5 Excitation waveform for 2D Cartesian and 3D rotationally symmetric solitons

One possible source/excitation waveform for 3D CRS solitons can be taken from 2D Cartesian solitons, where two parallel 2D solitons are $\mathrm{\pi }$ out of phase (antiphase) of each other. The 2D solitons are then modified to be cylindrical for the 3D CRS case, thereby generating a 3D wave as shown in Fig. 1.

For the case of a 3D CRS soliton, the 3D excitation is the same in the y = 0 plane as the two parallel antiphase soliton excitations from 2D. The two parallel antiphase soliton source from 2D is described in [5], and a modified version of the 2D source is shown below in Eq. (28) with a $\mathrm{\pi }$ phase difference between solitons, where the center of the excitation waveform is at x = 0,

A “source spacing” (spac) is defined by its usage in Eq. (28), where ${{x}_{{spac}}} = ({{w}_{0}}\; {spac) }$ for 2D and the distance between the two parallel antiphase soliton extremums is ${2}{{x}_{{spac}}}$, ${{w}_{0}}$ is the characteristic width, ${{w}_{c}}$ is the input carrier angular frequency and ${{A}_{0}}$ is amplitude of the input magnetic field.

The input waveforms ${{H}_{y}}({z = 0,y = 0} )$ plotted in Fig. 2 show excitation waveforms for a variety of “source spacings”. From the excitation waveforms shown in Fig. 2(a)), for source spacings roughly below spac = 4, the peak of ${{H}_{y}}{(\textrm{z}\; = \; 0,y\; = \; 0)}$ is seen to be less one, and must be corrected if ${{A}_{0}}$ is the desired peak value of the excitation waveform ${{H}_{y}}{(\textrm{z}\; = \; 0,\; y\; = \; 0)}$ as shown by Fig. 2(b), where the “${{A}_{0}}$ Correction Factor” is plotted for a number of source spacings.

For all of the simulations in this paper, ${{A}_{0}}$ has been corrected ($[{{{A}_{0}}{\; {Correction}\; {Factor}}} ]\mathrm{\ast }{{A}_{0}}$) as shown in Fig. 2(b) so that all source maximum values are the same regardless of source spacing as shown in Fig. 3, allowing for the easiest comparison of results of different source spacings.

 

Fig. 2. The excitation waveform ${{H}_y}{(z\; = \; 0,y\; = \; 0)}$ with a fixed ${{A}_{0}}$ and ${\sin}({{{w}_{c}}{t}} ){ = 1}$. (a) Shows the excitation waveform for different source spacings and the drop in the amplitude of the waveform maximum with smaller source spacings. (b) The ${{A}_{0}}$ Correction Factor required to keep the maximum ${{H}_y}{(z\; = \; 0,y\; = \; 0)\; = \; 1}$ regardless of source spacing.

Download Full Size | PDF

While the 2D source is only valid in the y = 0 plane ($\phi \mathrm{\;\ =\ \;\ 0^\circ }$ and $\phi \mathrm{\;\ =\ \;\ 180^\circ }$) at z = 0, the expanded 3D CRS source, obtained by observation, combining Eq. (28) and the rotational symmetry requirements, is valid for all angles of $\phi $ at z = 0. The 3D CRS source equation is defined below in Eq. (29), and the resulting ${{H}_y}{(z\; = \; 0)}$ in the y = 0 plane from Eq. (29) is plotted in Fig. 3, where ${{H}_\phi }{(z\; = \; 0,\; \textrm{r}\; \; } \ge {\; 0)}$ is shown at $\phi \mathrm{\;\ =\ \;\ 0^\circ }$ and $- {{H}_\phi }{(z\; = \; 0,\; \textrm{r}\; \; } \ge {\; 0)}$ is shown at $\phi \mathrm{\;\ =\ \;\ 180^\circ }$. A “source spacing” (spac) in this case, is defined by its usage in Eq. (29) below, where ${{r}_{{spac}}}{ = \; (}{{w}_{0}}{\; {spac})}$ for 3D and the distance between the hyperbolic secant extremums is ${2}{{r}_{{spac}}}$, and ${r\; = \; }\sqrt {{{x}^{2}}{ + }{{y}^{2}}} $ (when $\phi {\; = \; 0}$ then ${y\; = \; 0}$. causing ${{r}_{{spac}}}{ = \; }{{x}_{{spac}}}$),

The 3D rotationally symmetric aspect of the source in Eq. (29) is shown in Fig. 1 and in Fig. 3 where the ${{H}_\phi }({z = 0,r \ge 0} )$ source equation is the same for all values of $\phi $. ${{H}_\phi }({z = 0,r \ge 0} )$ plotted in the y = 0 plane, highlights the changing direction of the cylindrical unit vector $\hat{\phi }$ with angle $\phi $, where ${{H}_\phi }({z = 0,r \ge 0,\; \phi \mathrm{\;\ =\ \;\ 180^\circ }} )$ end up pointing in the opposite direction of ${{H}_\phi }({z = 0,r \ge 0,\; \phi \mathrm{\;\ =\ \;\ 0^\circ }} )$. The ${{H}_\phi }({z = 0,r \ge 0,\; \phi \mathrm{\;\ =\ \;\ 0^\circ }} )$ and ${{H}_\phi }({z = 0,r \ge 0,\; \phi \mathrm{\;\ =\ \;\ 180^\circ }} )$ fields together are equivalent to the 2D source ${{H}_{y}}{(z\; = \; 0)}$ in the special case.

 

Fig. 3. ${{H}_{y}}{(z\; = \; 0)}$ versus the normalized radius in the y = 0 plane resulting from the CRS source ${{H}_\phi }{(z\; = \; 0,r\; \; } \ge {\; 0,\; }\phi {)}$ with a “Corrected” ${{A}_{0}}$ and ${\sin}({{{w}_{c}}{t}} ){ = 1}$ for different source spacings.

Download Full Size | PDF

3. Numerical results and analysis

A systematic study was performed, simulating solitons for a series of “source spacings” from 1 to 15 with a set maximum value of the excitation waveform, using the material parameters and a selection of simulation parameters from [2,3,5] for the soliton configuration. Both 3D CRS and 2D Cartesian solitons were simulated for each “source spacing” from 1 to 5, to compare the results, and attempt to separate the 3D behavior from the 2D behavior.

All the computational simulation results in this article were conducted with the advanced computing resources provided by Texas A&M High Performance Research Computing.

3.1 Material and simulation parameters

The material parameters used in the simulations below were taken from [2,3] and are shown in Table 1. The input amplitude (${{A}_{0}}$) corrected for source spacing, characteristic width (${{w}_{0}}$) and input carrier frequency (${{w}_{c}}$) required for creating a 2D soliton were used as well, as shown in Table 2 along with other simulation parameters.

Table 1. Material Parameters^a

View Table | View all tables in this article

Table 2. Simulation Parameters^a

View Table | View all tables in this article

The 3D $\textrm{T}{\textrm{M}_\textrm{z}}$ CRS fields (${\textrm{H}_\mathrm{\phi} }$, ${\textrm{E}_\textrm{r}}$, ${\textrm{E}_\textrm{z}}$) were simulated using a nr x nz dimension Lebedev grid [12–14]. Liu’s paper [14] calls this grid the “unstaggered grid”, which is a combination of two shifted Yee grids with collocated electric field values. The radial spatial increment $\mathrm{\Delta} r$ and the total radial distance were chosen to minimize reflections from the boundary. A 3^rd order Newton backward difference polynomial interpolation was used for ${\textrm{H}_\mathrm{\phi} }$ as a basic boundary condition and a further safeguard against reflections. The soliton wavefront was never allowed to reach the end of the simulation space to avoid any reflections from the horizontal boundary. The simulation parameters $\mathrm{\Delta} z$, $\mathrm{\Delta} r$ and $\mathrm{\Delta} t$ were chosen empirically by performing numerical studies, and the resulting time-step is 7.527 times the required time-step for the Courant stability limit. For 2D $\textrm{T}{\textrm{M}_\textrm{z}}$ simulations, source symmetry, even symmetry for ${\textrm{E}_\textrm{z}}$ and odd symmetry with ${\textrm{H}_\mathrm{\phi} }$ (${\textrm{H}_\textrm{y}}$) and ${\textrm{E}_\textrm{r}}$ (${\textrm{E}_\textrm{x}}$), were utilized as boundary conditions at x = r = 0, allowing the same nr x nz dimension Lebedev grid to be used.

3.2 Simulation results and analysis

The simulation results for source spacing = 1 and source spacing = 5 are shown for both the 2D and the 3D CRS cases in Fig. 4 in the y = 0 plane. These source spacings were chosen to show the differences between 2D and 3D. It is interesting that both 2D and 3D simulations with source spacing = 1, show the soliton repelling/expanding, moving away from r = 0, as the wave propagates in the z direction. For source spacing = 5, the 2D solitons appear to propagate without attracting or repelling, while the 3D CRS soliton collapses towards r = 0.

 

Fig. 4. Results of |E| for 2D and 3D CRS soliton simulation for source spacings of 1 and 5. (a) 3D CRS with source spacing = 1, (b) 2D with source spacing = 1, (c) 3D CRS with source spacing = 5 and (d) 2D with source spacing = 1. Note: Fig. 4(c) uses the same simulation results as Fig. 1, where the 3D aspect of the CRS waveforms is shown.

Download Full Size | PDF

To compare the 3D CRS and 2D soliton propagation in the y = 0 plane, curves are constructed as shown in Fig. 5, starting at the ${{H}_\phi }$ extremum location in the r direction nearest the source location z = 0, then moving in the + z-direction connecting the subsequent ${{H}_\phi }$ soliton sinusoidal extremum locations in the r direction, where the extremum locations in the z direction result from the source equation term ${\sin(}{{w}_{c}}\mathrm{t)\ =\ \;\ \pm 1}$. The extremum locations are plotted versus propagation distance z in Fig. 5(b) for a series of source spacings from 1 to 5. By comparing the 2D and 3D solitons from Fig. 4 with Fig. 5(b) it can be seen that the soliton sinusoidal extremum locations of ${{H}_\phi }$ in Fig. 5(b) for source spacing = 1 and 5, match the soliton extremums in Fig. 4. It is interesting to note that the 3D CRS soliton expands like the 2D solitons repel for source spacings smaller than 2, appearing to keep that aspect of 2D behavior, and then collapses for spacings larger than 2. For 2D, with source spacings greater than 3, it appears the soliton repelling force is negligible; while for 3D CRS, with source spacings greater than 2, the soliton collapses towards r = 0. Soliton collapse is a distinctly 3D CRS behavior not observed in 2D.

 

Fig. 5. Soliton extremum locations of ${{H}_\phi }$ in the r direction versus propagation direction z in the y = 0 plane. (a) 3D CRS for source spacings from 1 to 15. (b) For comparison, 2D and 3D CRS for source spacings from 1 to 5. Due to the large extremum location density, the plot markers are at evenly sampled locations, while maintaining the shape of the curves.

Download Full Size | PDF

Next, the displacement in the r direction and slope of the soliton extremums from their starting locations nearest z = 0 is determined,

The displacement in r of the soliton sinusoidal extremum locations of ${{H}_\phi }$ in the transverse direction r from the starting soliton sinusoidal extremum locations of ${{H}_\phi }$ in the transverse direction r verses propagation distance z, is shown in Fig. 6(a). For each of the displacement curves shown in Fig. 6(a), the average slope is determined, and then plotted versus source spacing for 2D, 3D CRS and the 3D CRS slope minus the 2D slope, in order to isolate the 3D CRS behavior from the 2D repelling/expanding behavior, as shown in Fig. 6(b). The average slope is determined by averaging the slope for each extremum location on the displacement curve,

A few observations can be made from Fig. 6 related to the slope of the displacement curves. First, all the 2D displacement curve slopes appear to be positive, as expected for two antiphase 2D Cartesian solitons which are known to repel each other [6,15] producing a positive slope, while the 3D CRS displacement curves include positive and negative slopes. The reason why negative slopes are present for some of the 3D CRS results, corresponding to the wave collapsing towards r = 0, is explained in Section 3.3. Second, for source spacings less than 2, the 3D CRS displacement curve slopes are similar in shape to the 2D displacement curves. The smaller the source spacing, the closer the 3D CRS displacement curve slope is to the 2D displacement curve slope, implying the 3D mimics the 2D repelling/expanding behavior for smaller source spacings. Third, 3D CRS source spacing = 2 seems to be near a transition point from positive average slope to negative average slope. Fourth, for 3D CRS, the max negative slope increases from being near zero at source spacing = 2 to reaching a maximum negative slope at source spacing = 3. Fifth, for 2D, the average slope approaches zero as the source spacing increases beyond source spacing = 3, where the average slope is effectively zero for source spacing = 5. Lastly, in an attempt to isolate the 3D specific collapsing behavior, looking at “3D” and “3D-2D” for larger than source spacing = 3, it is seen that the max negative slope decreases in magnitude in a roughly -1/(spac) shape, decreasing towards zero as source spacing increases. This shows that as the source spacing increases to larger than 5, the 3D non-zero slope is approaching the known zero slope behavior of two largely separated antiphase 2D solitons, as predicted in the discussion around Eq. (24).

 

Fig. 6. (a) The displacement in r from the ${{\textrm{H}}_\phi }$ extremum location nearest the excitation at z = 0 to the ${{\textrm{H}}_\phi }$ extremum location in r at each evenly sampled z extremum location, shown in Fig. 5. (b) The “Average Slope”, as can be seen visually from Fig. 6(a), for 2D Average Slope, 3D Average Slope and [3D Average Slope - 2D Average Slope] vs. source spacings.

Download Full Size | PDF

3.3 Interpretation of simulation results behavior

From the simulation results, it appears the longitudinal fields and polarization currents play a key role in the cause and effect of the collapsing and expanding behavior/process. This can be seen by the overall picture presented by comparing the 3D CRS results of this section with the 2D results of this section.

An argument for “why” the 3D CRS wave collapses is related to the “signs” of the two terms in the equation ${({\nabla {\ \times \textbf{H}}} )_{{\textrm{z,3D}}}}{ = }\frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}$, since ${({\nabla \mathrm{\ \times \textbf{H}}} )_{{\textrm{z,3D}}}}$ is one of the dominant terms used in determining ${{E}_{z}}$ using the Newton-Raphson method, which is then used, along with ${{E}_{r}}$, to determine ${{H}_\phi }$ using the update equation. The source waveform for source spacing = 5 is shown in Fig. 7 with labels for the “Inner Region” and the “Outer Region” for both ${\sin}({{{w}_{c}}{t}} ){ = 1}$ and ${\sin}({{{w}_{c}}{t}} ){ = } - {1}$.

 

Fig. 7. Slope of ${{H}_\phi }$ ($\partial {{\textrm{H}}_\mathrm{\phi} }/\partial {\textrm{r}}$) and ${{H}_\phi }/r$ sign comparison plot corresponding to the equation ${({\nabla {\textrm{xH}}} )_{{\textrm{z,3D}}}}{ = }\frac{{\partial {{\textrm{H}}_\mathrm{\phi} }}}{{\partial {\textrm{r}}}}{ + }\frac{{{{\textrm{H}}_\mathrm{\phi} }}}{{\textrm{r}}}$ for source spacing = 5. (a) ${\sin}({{{w}_{c}}{t}} ){ = 1}$. (b) ${\sin}({{{w}_{c}}{t}} ){ = } - {1}$.

Download Full Size | PDF

In the “Inner Region”, ${{H}_\phi }$ is the same sign as the slope of ${{H}_\phi }$ (either both positive or both negative); this causes ${({\nabla {\ \times \textbf{H}}} )_{{\textrm{z,3D}}}}$ to grow in the “Inner Region”, effectively steepening the inner slope of the ${{H}_\phi }$ as it propagates. In the “Outer Region”, ${{H}_\phi }$ is the opposite sign as the slope of ${{H}_\phi }$; this causes ${({\nabla {\ \times \textbf{H}}} )_{{\textrm{z,3D}}}}$ to decrease in the “Outer Region”, effectively flattening the outer slope of the ${{H}_\phi }$ as it propagates.

For the 3D CRS case, with source spacing = 5, the cross section simulation results for ${{H}_\phi }$, $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}$ and $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}$ are plotted at a selected z location as shown in Fig. 8 to emphasize the difference between 2D and 3D CRS cases, as an example of the results of the effect discussed in Fig. 7.

 

Fig. 8. ${{H}_y} = {{H}_\phi }$ at $\phi \mathrm{\;\ =\ \;\ 0^\circ }$, ${{H}_y} ={-} {{H}_\phi }$ at $\phi \mathrm{\;\ =\ \;\ 180^\circ }$, $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}$ and $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}$ at the ${z\; = \; 8}\mathrm{.112\;\ \mu m}$ cross section for the 3D CRS simulation with source spacing = 5. This shows the additive behavior in the “Inner Region”, and the subtractive behavior in the “Outer Region” for the ${({\nabla \mathrm{\ \times \textbf{H}}} )_{z,3D}} = \frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}$ term as discussed in Fig. 7 compared to the $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}$ term which corresponds to the 2D term in the y = 0 plane, where the $\frac{{\partial {{H}_\phi }}}{{\partial {r}}}$ curve is also the ${({\nabla \mathrm{\ \times \textbf{H}}} )_{z,2D}} = \frac{{\partial {{H}_y}}}{{\partial {x}}}$ curve for the plotted ${{H}_y}$.

Download Full Size | PDF

At the selected z location, the cross-section of the terms used in the Newton-Raphson method are plotted as shown in Fig. 9 and Fig. 10 for the 3D CRS and 2D cases. The effect described above can be seen in Fig. 8, where $\left|{\frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}} \right|$ is greater than $\left|{\frac{{{{H}_\phi }}}{{r}}} \right|$ in the “Inner Region”, and $\left|{\frac{{\partial {{H}_\phi }}}{{\partial {r}}}{ + }\frac{{{{H}_\phi }}}{{r}}} \right|$ is less then $\left|{\frac{{{{H}_\phi }}}{{r}}} \right|$ in the “Outer Region”. The effect can also be seen in the curves shown in Fig. 9(a),(c), about the maximum of ${{E}_{r}}$ and the zero crossing of ${({\nabla {\ \times \textbf{H}}} )_{z}}$, where the Newton-Raphson terms for the 3D CRS case have larger magnitude values in the “Inner Region” than in the “Outer Region”. For 2D, it can be seen in Fig. 10(a),(c), about the maximum of ${{E}_{r}}$ and the zero crossing of ${({\nabla {\ \times \textbf{H}}} )_{z}}$, where the Newton-Raphson terms for the 2D case have roughly equal magnitude values in the “Inner Region” and in the “Outer Region”. This shows a “quasi symmetry” in the 2D case which is not present in the 3D CRS case.

 

Fig. 9. 3D CRS, ${{E}_{z}}$ and ${{E}_{r}}$ along with the corresponding Newton-Raphson method Lorentz, Kerr and Raman polarization currents, time derivative and the curl, from Eq. (7) at the ${z\; = \; 8}\mathrm{.112\;\ \mu m}$ cross section with source spacing = 5. ${({\nabla \mathrm{\ \times \textbf{H}}} )_{z}}$ equation terms on (a) a linear plot and (b) their magnitude on a semilog plot. ${({\nabla \mathrm{\ \times \textbf{H}}} )_r}$ equation terms on (c) a linear plot and (d) their magnitude on a semilog plot.

Download Full Size | PDF

 

Fig. 10. 2D, ${{E}_{z}}$ and ${{E}_{r}}$ along with the corresponding Newton-Raphson method Lorentz, Kerr and Raman polarization currents, time derivative and the curl, from Eq. (5) at the $z = 8.112\mathrm{\;\ \mu m}$ cross section with source spacing = 5. ${({\nabla \mathrm{\ \times \textbf{H}}} )_{z}}$ equation terms on (a) a linear plot and (b) their magnitude on a semilog plot. c-d) ${({\nabla \mathrm{\ \times \textbf{H}}} )_r}$ equation terms on (c) a linear plot and (d) their magnitude on a semilog plot.

Download Full Size | PDF

The moment of the magnitude of the sum of the Lorentz, Kerr, and Raman polarization current waveforms has been calculated using the clearly weighted and unsymmetric polarization current distributions shown in Fig. 9 and Fig. 10, about the ${{E}_r}$ extremum r location and ${({\nabla {\ \times \textbf{H}}} )_{z}}$ zero crossing r location. In mechanics a “moment” is the tendency of a body to rotate about an axis, given an applied distributed force. In our case, we will attempt to use the moment to quantify the tendency for the wave to collapse or expand (attract or repel in 2D) given a polarization current distribution in the r direction.

The moment of $|{\textrm{J}_{{r/z}}}{(\textrm{r},\; \textrm{z})|}\; $“about the ${\textrm{H}_\mathrm{\phi} }$ extremum location in r” for a given z location with ${\mathrm{\rho }_{{\textrm{H}_\mathrm{\phi} }{\; \textrm{extremum}\; \textrm{location}}}}$, is shown in Eq. (33) and Eq. (34),

For each of the selected z locations, the moments of the polarization current, $|{\textrm{J}}_{z}|= |{\textrm{J}}_{{z,\textrm{Lorentz}}}+{\textrm{J}}_{{{z,\textrm{Kerr}}}}+ {\textrm{J}}_{{{z,\textrm{Raman}}}}|$ and ${|}{{\textrm{J}}_{r}}{|= |}{{\textrm{J}}_{{{r,\textrm{Lorentz}}}}}{ + }{{\textrm{J}}_{{r,\textrm{Kerr}}}}{ + }{{\textrm{J}}_{{r,\textrm{Raman}}}}{|}$ waveforms have been calculated. All moments in this paper have been calculated as defined by Eq. (33) and Eq. (34), “about the ${{H}_\phi }$ extremum location in the r direction” plotted in Fig. 5. For example, at ${z\; = \; 8}\mathrm{.112\mu m}$, to find the moment of $|{{{\textrm{J}}_{z}}} |$, the ${{\textrm{J}}_{z}}$ shown in Fig. 9(a) would be used, and then to find the moment of $|{{{\textrm{J}}_r}} |$, the ${{\textrm{J}}_r}$ shown in Fig. 9(c) would be used. The moments of $|{{{\textrm{J}}_r}} |$ calculated at ${{H}_\phi } $ extremum z locations versus propagation distance z are shown in Fig. 11 for selected source spacings. The moments of $|{{{\textrm{J}}_z}} |$ calculated at ${{H}_\phi } $ extremum z locations versus propagation distance z are shown in Fig. 12(a) for selected source spacings. And lastly, the smoothed sum of the moment of ${|}{{\textrm{J}}_{z}}{|}$ and the moment of ${|}{{\textrm{J}}_{r}}{|}$ at ${{H}_\phi } $ extremum z locations versus propagation distance z are shown in Fig. 12(b) for the same selected source spacings.

 

Fig. 11. Normalized and smoothed moments for 4 spacings using a moving average for ${|}{{J}_r}{|}$ about ${{H}_\phi }$ along the z-axis at extremum locations in r, emphasizing the oscillatory nature of the ${|}{{J}_r}{|}$ moment while showing the overall trend. (a) 3D CRS simulations and (b) 2D simulations.

Download Full Size | PDF

 

Fig. 12. (a) Normalized Moment of ${|}{{J}_{z}}{|}$ about ${{H}_\phi }$ along the z-axis at evenly sampled extremum location in r for 2D and 3D CRS simulations. (b) Normalized Smoothed([Moment of ${|}{{J}_{r}}{|}$ about ${{H}_\phi }$ extremum location in r] + [Moment of ${|}{{J}_{z}}{|}$ about ${{H}_\phi }$ extremum location in r]) for 2D and 3D CRS simulations.

Download Full Size | PDF

A clear difference between the 2D and 3D CRS moments of the longitudinal polarization currents ${|}{{J}_{z}}{|}$ can be seen in Fig. 12(a). The 3D CRS moments of ${|}{{J}_{z}}{|}$ are mostly negative, while all the 2D moments of ${|}{{J}_{z}}{|}$ are always positive. The negative 3D or positive 2D moment of ${|}{{J}_{z}}{|}$ behavior generally correlates with the observed soliton behavior in Fig. 5: for 3D, a tendency towards collapsing, with moderate source spacings collapsing the fastest and the rate of collapsing decreasing for larger source spacings as shown by source spacing = 15; and for 2D, a tendency towards repelling for smaller source spacings, with no interaction between the antiphase solitons for larger source spacings. It should be noted that the moment of ${|}{{J}_{z}}{|}$ near the front of the soliton, z > 25$\mathrm{\mu}$m, does not seem to correlate as closely with the behavior of the wave propagation, as can be seen by comparing with Fig. 5. That said, for z < 25$\mathrm{\mu}$m, it appears the ${{J}_{z}}$ distribution in the r direction strongly correlates with the behavior of the wave propagation for 2D and 3D, and may slowly impact it. There is also a clear correlation between the wave propagation and the moments of the transverse polarization currents ${|}{{J}_{r}}{|}$ shown in Fig. 11 along with the sum of the moments of the transverse and longitudinal polarization currents shown in Fig. 12(b).

In the optical spatial soliton literature review article by Stegeman and Segev [15], the effective refractive index value is said to be “intensity dependent” and the cause of two parallel solitons repelling or attracting depending on the phase. As the wave propagates, the intensity dependent effective index of refraction decreases or increases in some regions, as the intensity of the wave decreases or increases, causing the wave to move away from smaller effective refractive index values and towards the larger effective refractive index values. Our results indicate that the polarization current distribution, shown in Fig. 9 and Fig. 10 as examples, in the r direction potentially guides the wave in a similar manner, moving towards the larger polarization current values. In the case of 3D CRS solitons, the extra term in the ${({\nabla {\ \times \textbf{H}}} )_{z}}$ causes larger magnitude polarization current values to exist in the “Inner Region” than in the 2D case, making the ${({\nabla {\ \times \textbf{H}}} )_{z}}$ term play a more prominent role. These larger longitudinal polarization currents created by ${({\nabla {\ \times \textbf{H}}} )_{z}}$ in the “Inner Region” correlate with and potentially impact the new collapsing behavior, either re-enforcing or countering the collapsing behavior.

4. Conclusion

In this paper we have presented, to our knowledge for the first time, simulation of 3D full wave time domain electromagnetic spatial soliton propagating in a nonlinear dispersive media with radial polarization. A new method for construction and analysis of a $\textrm{T}{\textrm{M}_\textrm{z}}$ rotationally symmetric 3D electromagnetic wave was developed using a two-dimensional numerical grid using the FDTD GVADE method modified for the cylindrical rotationally symmetric problem. The behavior of the 3D hyperbolic secant cylinder soliton, composed of two subtracted radial hyperbolic secant functions, was compared with two parallel antiphase 2D hyperbolic secant solitons and therefore physically similar in cross-section to the 3D cylinder, allowing analysis to be performed and conclusions to be drawn between these two dimensionally different but visually similar soliton waves. An important result is the tendency of the 3D cylinder to collapse towards r = 0, with moderate source spacings collapsing the fastest and the rate of collapsing decreasing as approximately 1/(source spacing) for larger source spacings. A possible explanation for why the 3D rotationally symmetric soliton collapses towards r = 0, was given. The theory and results of this paper allow for a possible improved understanding of some 2D behavior, using polarization current distributions in the Newton-Raphson method, and insights into the behavior of 3D cylindrical rotationally symmetric waves in a nonlinear medium.