PDF

Self-Focusing
Introduction
For low intensities, the refractive index is independent of the intensity of the light in the material. However, when the intensity is large, so large that the electric field of the light field actually start to perturb the electron clouds around the nuclei, the refractive index start to depend on the intensity. Thus, for dielectrics, like glass, the refractive index increases with the intensity.
When a Gaussian beam propagates through a medium with an intensity-dependent refractive index, the index is highest at the center of the beam. Thus, the induced refractive index profile acts as a lens or a waveguide that counteracts the spreading of the beam due to diffraction. The effect that the beam itself induces this positive, focusing, lens in the material is called self-focusing.
Self-focusing manifests itself both in terms of whole-beam focusing, where the beam’s properties are changed by the induced index profile, and by small-scale focusing, where noise across the beam’s cross-sectional intensity distribution is amplified and the beam can break up in to several self-focused filaments.
The nonlinear refractive index is written as
where n0 is the constant (linear) part of the refractive index, γ is the nonlinear refractive index coefficient and I is the intensity. The nonlinearity is due to the optical Kerr effect, which is a nonlinear distortion of the electron clouds around the nuclei in the material. For the standard optical glass BK-7, the nonlinear coefficient is 4·10-16 cm2/W. Thus, for an intensity of 2.5 GW/cm2, the induced refractive index change is 10-6. This is a small change, but it has a significant effect on the beam parameters.
Self-focusing is important from a laser engineering perspective, as the modification of the beam must be incorporated in the design. Furthermore, if the threshold for self-focusing is exceeded, the material is damaged. Thus, it is important to know the self-focusing threshold values for the materials used in the design. Self-focusing occurs in dielectrics, like optical glasses and laser rod materials, such as Nd:YAG.
A first estimate of the threshold power for self-focusing is obtain by assuming that the beam has a circular cross-section with a uniform intensity. Within this beam, the refractive index is higher than outside the beam. Thus, the beam itself induces a waveguide structure. You can equal the acceptance angle of the waveguide with the diffraction angle from the circularly confined light. From this equality, you get the critical power to be
(1)
Model Definition
The geometry for the model is simple - just a cylinder.
The incident beam is approximated by a Gaussian beam1, polarized in the z-direction. The electric field is given by
where E0 is the electric field amplitude, w0 is the spot radius at the waist, k is the wave number, defined by k = 2πn0/λ, and z is the unit vector in the z-direction. The function w(x) defines the spot size variation as a function of the distance from the beam waist,
(2)
where x0 = πn0w02/λ is the Rayleigh range. The radius of curvature is defined by
and the phase change close to the beam waist, the so-called Gouy shift, is defined by
Results and Discussion
Figure 1 shows the Gaussian beam for a low input intensity. As expected, the beam is symmetric around the beam waist location.
Figure 1: The Gaussian beam for a low peak intensity, I0 = 1 kW/cm2.
Figure 2 shows the beam for a high input intensity, I0 = 14 GW/cm2. The nominal induced refractive index, γI0, is 5.6×10-6. As noted from Figure 3 the actual induced refractive index is much larger than the nominal value. This is of course an effect of the self-focusing of the beam.
Figure 2: The Gaussian beam for a high peak intensity, I0 = 14 GW/cm2.
To demonstrate the effect of self-focusing, Figure 4 shows the calculated spot radius at the boundary between the propagation domain and the PML domain. The spot radius is defined by
(3)
where A is the integration area and I(y,z) is the cross-sectional intensity distribution of the beam. For low peak intensities, the calculated spot radius, using Equation 3, give similar results as the Gaussian beam expression in Equation 2. However, with increasing intensities, the beam start to deviate from a Gaussian beam and the spot radius is reduced linearly with intensity.
Figure 3: The induced refractive index change, γI, for a high peak intensity, I0 = 14 GW/cm2.
Figure 4: The spot radius at the end of the propagation domain versus the peak intensity.
The final intensity used in the parametric sweep is 14 GW/cm2. This corresponds to a power that is 23% of the critical power, provided in Equation 1. As discussed in Ref. 1, it is expected that the critical power for a Gaussian beam is reduced by a factor of approximately 4/(1.22π)2 = 0.27. Thus, the power corresponding to the last peak intensity in the sweep is approximately a factor 0.23/0.27 = 0.85 from the critical power for self-focusing for a Gaussian beam, where the induced refractive index profile completely balances the diffractive spreading of the beam.
To compute the field for even higher intensities, a much finer mesh would be needed, as the beam can break up into filaments. For a nonlinear problem like this, it is important to verify the results by, for instance, repeating the simulation with a finer mesh.
Reference
1. W. Koechner, Solid-State Laser Engineering, Springer, chap. 4.6, 2010.
Application Library path: Wave_Optics_Module/Nonlinear_Optics/self_focusing
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Optics>Wave Optics>Electromagnetic Waves, Beam Envelopes (ewbe).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Wavelength Domain.
6
Global Definitions
Start by adding some global model parameters.
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Geometry 1
Cylinder 1 (cyl1)
1
In the Geometry toolbar, click  Cylinder.
2
In the Settings window for Cylinder, locate the Size and Shape section.
3
In the Radius text field, type radius.
4
In the Height text field, type length.
5
Locate the Position section. In the x text field, type -length/2.
6
Locate the Axis section. From the Axis type list, choose x-axis.
Block 1 (blk1)
Since the planes y = 0 and z = 0 should be symmetry planes for the beam, we only need to model one forth of the beam. Thus, add a block that intersect one quadrant of the cylinder.
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type length.
4
In the Depth text field, type radius.
5
In the Height text field, type radius.
6
Locate the Position section. In the x text field, type -length/2.
7
In the y text field, type -radius.
Intersection 1 (int1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Intersection.
2
Click in the Graphics window and then press Ctrl+A to select both objects.
3
In the Geometry toolbar, click  Build All.
Global Definitions
Since the geometry is so long and thin, modify the view settings to not preserve the aspect ratio.
Definitions
In the Model Builder window, expand the Component 1 (comp1)>Definitions node.
Camera
1
In the Model Builder window, expand the Component 1 (comp1)>Definitions>View 1 node, then click Camera.
2
In the Settings window for Camera, locate the Camera section.
3
From the View scale list, choose Automatic.
4
From the Automatic list, choose Anisotropic.
5
Select the Automatic update check box.
6
Click  Update.
7
Click the  Zoom Extents button in the Graphics toolbar.
Materials
Now add the BK-7 glass used in the model.
BK-7 glass
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, type BK-7 glass in the Label text field.
3
Locate the Material Contents section. In the table, enter the following settings:
The variable ewbe.Poavx represents the intensity in the propagation direction.
Definitions
Setup a boundary integration operator, for calculation of the output power and the output spot radius.
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, type intop_output_boundary in the Operator name text field.
3
Locate the Source Selection section. From the Geometric entity level list, choose Boundary.
4
Variables 1
Add the expressions for the output power and the output spot radius.
1
In the Definitions toolbar, click  Local Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Electromagnetic Waves, Beam Envelopes (ewbe)
Set the interface to use unidirectional propagation and define the wave vector component in the x direction.
1
In the Model Builder window, under Component 1 (comp1) click Electromagnetic Waves, Beam Envelopes (ewbe).
2
In the Settings window for Electromagnetic Waves, Beam Envelopes, locate the Wave Vectors section.
3
From the Number of directions list, choose Unidirectional.
4
Specify the k1 vector as
Perfect Magnetic Conductor 1
Add a perfect magnetic conductor boundary condition on the symmetry plane y = 0. Since the other symmetry plane, z = 0, is an exterior boundary, there is already a default perfect electric conductor boundary condition applied on that boundary.
1
In the Physics toolbar, click  Boundaries and choose Perfect Magnetic Conductor.
2
Use a matched boundary condition to launch an incident Gaussian beam polarized in the z direction.
Matched Boundary Condition 1
1
In the Physics toolbar, click  Boundaries and choose Matched Boundary Condition.
2
3
In the Settings window for Matched Boundary Condition, locate the Matched Boundary Condition section.
4
From the Incident field list, choose Gaussian beam.
5
In the w0 text field, type w0.
6
In the p0 text field, type length/2.
7
From the Input quantity list, choose Power.
8
In the P text field, type pi*w0^2/2*I0.
9
Specify the Eg0 vector as
Add a reference point that together with the propagation direction defines the optical axis.
Reference Point 1
1
In the Physics toolbar, click  Attributes and choose Reference Point.
2
Select Point 2 only, that is the point on the optical axis (where y = 0 and z = 0). The focal point is located the distance p0 from this reference point in the propagation direction, defined by ewbe.k1.
Matched Boundary Condition 2
1
In the Physics toolbar, click  Boundaries and choose Matched Boundary Condition.
2
Mesh 1
The physics-controlled meshing cannot handle the nonlinear refractive index defined for the BK-7 material. Thus, create a custom mesh. First create a triangular mesh on the input boundary and then a swept mesh for the domain.
Free Triangular 1
In the Mesh toolbar, click  Boundary and choose Free Triangular.
Size
Set the main Size settings to represent the discretization along the propagation direction. Since the beam is expected to behave as a slightly distorted Gaussian beam, it is sufficient to set the maximum mesh element size to half the Rayleigh range.
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type x0/2.
Free Triangular 1
1
In the Model Builder window, click Free Triangular 1.
2
Size 1
1
Right-click Free Triangular 1 and choose Size.
The triangular mesh should resolve the beam’s cross-sectional distribution.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. Select the Maximum element size check box.
5
6
Select the Minimum element size check box.
7
Create a swept mesh in the propagation direction.
Swept 1
1
In the Mesh toolbar, click  Swept.
2
In the Model Builder window, right-click Mesh 1 and choose Build All.
Study 1
Do not generate the default plots.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, locate the Study Settings section.
3
Clear the Generate default plots check box.
Step 1: Wavelength Domain
1
In the Model Builder window, under Study 1 click Step 1: Wavelength Domain.
2
In the Settings window for Wavelength Domain, locate the Study Settings section.
3
In the Wavelengths text field, type lda0.
Parametric Sweep
Setup a parametric sweep of the nominal peak intensity, from 1 kW/cm2 (corresponding to linear propagation) to 14 GW/cm2 that will show a significant self-focusing effect.
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
4
5
Click to expand the Advanced Settings section. From the Use parametric solver list, choose Off, to turn off the parametric solver that otherwise would perform calculations also for intermediate intensities.
Solution 1 (sol1)
Since this is a nonlinear problem, it is better to split the complex electric field variable into its real and imaginary parts. This will produce a more accurate Jacobian for the problem, leading to a faster convergence.
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node, then click Compile Equations: Wavelength Domain.
3
In the Settings window for Compile Equations, locate the Study and Step section.
4
Select the Split complex variables in real and imaginary parts check box.
5
In the Study toolbar, click  Compute.
Results
3D Plot Group 1
1
In the Home toolbar, click  Add Plot Group and choose 3D Plot Group.
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 1/Parametric Solutions 1 (sol2).
Slice 1
1
Right-click 3D Plot Group 1 and choose Slice.
2
In the Settings window for Slice, locate the Plane Data section.
3
From the Plane list, choose XY-planes.
4
From the Entry method list, choose Coordinates.
Deformation 1
1
Right-click Slice 1 and choose Deformation.
2
In the Settings window for Deformation, locate the Expression section.
3
In the Z component text field, type ewbe.normE.
3D Plot Group 1
1
In the Model Builder window, click 3D Plot Group 1.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges check box.
4
Click the  Go to Default View button in the Graphics toolbar.
Take a look at the field distributions for the different intensities.
5
Locate the Data section. From the Parameter value (I0 (W/m^2)) list, choose 1E7.
6
In the 3D Plot Group 1 toolbar, click  Plot.
Slice 1
Click the  Zoom Extents button in the Graphics toolbar. Compare the graph with that in Figure 1.
3D Plot Group 1
1
In the Model Builder window, click 3D Plot Group 1.
2
In the Settings window for 3D Plot Group, locate the Data section.
3
From the Parameter value (I0 (W/m^2)) list, choose 2E13.
4
In the 3D Plot Group 1 toolbar, click  Plot.
5
From the Parameter value (I0 (W/m^2)) list, choose 5E13.
6
In the 3D Plot Group 1 toolbar, click  Plot.
7
From the Parameter value (I0 (W/m^2)) list, choose 8E13.
8
In the 3D Plot Group 1 toolbar, click  Plot.
9
From the Parameter value (I0 (W/m^2)) list, choose 1.1E14.
10
In the 3D Plot Group 1 toolbar, click  Plot.
11
From the Parameter value (I0 (W/m^2)) list, choose 1.4E14.
12
In the 3D Plot Group 1 toolbar, click  Plot.
13
Click the  Zoom Extents button in the Graphics toolbar. Compare the graph with that in Figure 2.
Slice 1
To really see how the beam changes with intensity, you can also visualize this by running the 3D plots in the Player.
The first step would be to normalize the field solution with the nominal peak electric field.
1
In the Model Builder window, click Slice 1.
2
In the Settings window for Slice, locate the Expression section.
3
In the Expression text field, type ewbe.normE/E0.
Deformation 1
1
In the Model Builder window, click Deformation 1.
2
In the Settings window for Deformation, locate the Expression section.
3
In the Z component text field, type ewbe.normE/E0.
Animation 1
1
In the Results toolbar, click  Animation and choose File.
2
In the Settings window for Animation, locate the Target section.
3
From the Target list, choose Player.
4
Locate the Animation Editing section. From the Loop over list, choose I0.
5
Locate the Playing section. In the Display each frame for text field, type 1.
6
Right-click Animation 1 and choose Play.
If you check the Repeat box in the Playing settings, you can have the Player repeat the sequence.
3D Plot Group 2
Now create a new plot of the refractive index change induced by the beam.
In the Results toolbar, click  3D Plot Group.
Slice 1
1
Right-click 3D Plot Group 2 and choose Slice.
2
In the Settings window for Slice, locate the Plane Data section.
3
From the Plane list, choose XY-planes.
4
From the Entry method list, choose Coordinates.
5
Locate the Expression section. In the Expression text field, type ewbe.nxx-n0.
Deformation 1
1
Right-click Slice 1 and choose Deformation.
2
In the Settings window for Deformation, locate the Expression section.
3
In the Z component text field, type ewbe.nxx-n0.
3D Plot Group 2
1
In the Model Builder window, click 3D Plot Group 2.
2
In the Settings window for 3D Plot Group, locate the Plot Settings section.
3
Clear the Plot dataset edges check box.
4
Click to expand the Title section. From the Title type list, choose Manual.
5
In the Title text area, type Nonlinear refractive index: ewbe.nxx-n0.
6
In the 3D Plot Group 2 toolbar, click  Plot.
7
Click the Go to View 1 button in the Graphics toolbar.
8
Click the  Go to Default View button in the Graphics toolbar. Compare your graph with Figure 3.
Global Evaluation 1
Now, visualize how the spot radius decreases with intensity, using a table and a table plot.
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Expressions section.
3
4
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol2).
5
Click  Evaluate.
Table
1
Go to the Table window.
2
Click Table Graph in the window toolbar.
Results
Table Graph 1
1
In the Model Builder window, under Results>1D Plot Group 3 click Table Graph 1.
2
In the Settings window for Table Graph, locate the Data section.
3
From the x-axis data list, choose I0 (W/m^2).
4
From the Plot columns list, choose Manual.
5
In the Columns list, select Spot radius on output boundary (m).
1D Plot Group 3
1
In the Model Builder window, click 1D Plot Group 3.
2
In the Settings window for 1D Plot Group, locate the Axis section.
3
Select the Manual axis limits check box.
4
In the y minimum text field, type 1.5e-4.
5
In the y maximum text field, type 2.5e-4.
6
In the 1D Plot Group 3 toolbar, click  Plot. Your result should look similar to that in Figure 4.
Global Evaluation 2
Now, compare the result with the output spot radius of the ideal Gaussian beam.
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Expressions section.
3
4
Click  Evaluate.
Table
1
Go to the Table window.
You should find that the spot radius for the low intensity cases are close to that of the nominal beam.
Results
Finally, compare the total power in the beam with the critical power for self-focusing, as defined in Figure 4.
Global Evaluation 3
1
In the Results toolbar, click  Global Evaluation.
2
In the Settings window for Global Evaluation, locate the Expressions section.
3
Here we multiplied the power ratio with 4 since we integrated the power over only a fourth of the beam.
4
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol2).
5
Click  Evaluate.
Table
1
Go to the Table window.
You should find that the power ratio is approximately 23%.
 

1
Notice that the Gaussian beam expression used here is not a solution to Maxwell’s equations, but to the approximate paraxial wave equation. Thus, it should not be used for beams with spot radii of the same size as or smaller than the wavelength.