FDTD method for laser absorption in metals for large scale problems

The FDTD method has been successfully used for many electromagnetic problems, but its application to laser material processing has been limited because even a several-millimeter domain requires a prohibitively large number of grids. In this article, we present a novel FDTD method for simulating large-scale laser beam absorption problems, especially for metals, by enlarging laser wavelength while maintaining the material’s reflection characteristics. For validation purposes, the proposed method has been tested with in-house FDTD codes to simulate p-, s-, and circularly polarized 1.06 μm irradiation on Fe and Sn targets, and the simulation results are in good agreement with theoretical predictions. ©2013 Optical Society of America OCIS codes: (050.1755) Computational electromagnetic methods; (350.3390) Laser materials processing; (260.3910) Metal optics. References and links 1. K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwells equations in isotropic media,” IEEE Trans. Antenn. Propag. 14(3), 302–307 (1966). 2. C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express 18(20), 21427–21448 (2010). 3. K. Kitamura, K. Sakai, and S. Noda, “Finite-difference time-domain (FDTD) analysis on the interaction between a metal block and a radially polarized focused beam,” Opt. Express 19(15), 13750–13756 (2011). 4. S. Buil, J. Laverdant, B. Berini, P. Maso, J. P. Hermier, and X. Quélin, “FDTD simulations of localization and enhancements on fractal plasmonics nanostructures,” Opt. Express 20(11), 11968–11975 (2012). 5. C. Lundgren, R. Lopez, J. Redwing, and K. Melde, “FDTD modeling of solar energy absorption in silicon branched nanowires,” Opt. Express 21(S3), A392–A400 (2013). 6. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005). 7. H. Ki and J. Mazumder, “Numerical simulation of femtosecond laser interaction with silicon,” J. Laser Appl. 17(2), 110–117 (2005). 8. H. Li and H. Ki, “Effect of ionization on femtosecond laser pulse interaction with silicon,” J. Appl. Phys. 100(10), 104907 (2006). 9. W. M. Steen and J. Mazumder, Laser Material Processing, 4th ed. (Springer-Verlag, 2010). 10. M. Born and E. Wolf, Principles of Optics, 7th ed. (Cambridge University, 1999).


Introduction
Since the FDTD algorithm was first developed by Yee in 1966 [1], it has been extensively used for a variety of electromagnetic problems [2][3][4][5][6].Recently, there have been some efforts in the laser material processing community to use the FDTD method for simulating laser material interaction problems [7,8].Because accurate predictions on laser absorption in materials is by far the most important in better understanding of the processes, the direct solving of the Maxwell equations by the FDTD method seems like an ideal method.
However, due to its very stringent requirement on wavelength-dependent grid density, this method has been considered inappropriate for simulating light interaction with materials, because the domain size is extremely large compared to the wavelength of the light.In fact, Ki et al. [7,8] simulated femtosecond laser interaction with silicon targets using the FDTD method, but in order to reduce the total grid number, they employed the body-of-revolution (BOR) FDTD method assuming a radially polarized laser beam.Besides, the domain size was only about 50 μm in their studies.If a full three dimensional code needs to be applied to a typical laser manufacturing process, however, a conventional FDTD method is still far from a viable option.In a typical laser welding problem, for example, an Nd:YAG or CO 2 laser is generally irradiated on a metal plate that is at least several millimeters thick.In the case of a Nd:YAG laser, the typical wavelength is 1.06 μm, and assuming that 10 grids per wavelength is required and the plate thickness is 1 mm, the grid number in one dimension is roughly 10,000, which leads to ~10 12 grids in three dimensions.If a CO 2 laser beam is used, which has a wavelength of 10.6 μm, the total grid number can be decreased to ~10 9 , but even this grid number can be handled only by the most powerful supercomputing systems.
In this article, we propose a FDTD based scheme of simulating laser absorption in metals that can be used for relatively large scale problems, such as laser manufacturing processes.In this method, instead of using the original wavelength of the laser beam, a scheme to use a much enlarged wavelength has been developed, where the angle-dependent absorption characteristic is not altered or minimally changed.For validation purposes, we have tested the method with both the original Yee algorithm (or standard FDTD, hereafter) and an FDTD algorithm for dispersive media with the Drude model [6] (or dispersive FDTD, hereafter).We have also proposed a scheme that enables the use of the standard FDTD method for dispersive media by obtaining a new set of refractive index and extinction coefficient.Numerical tests have been performed for 1.06 μm laser beam interaction with iron (Fe) and tin (Sn) targets.The obtained simulation results are in good agreement with the theoretical predictions.

Changing wavelength for standard FDTD algorithm
Light absorption in a metal can be determined by the metal's complex refractive index and the incident angle as shown by the reflectance (R) formulas for s-and p-polarized lights where subscripts s and p denote s-and p-polarizations, n and κ are the real and imaginary parts of the complex refractive index n  (i.e., refractive index and extinction coefficient), and i θ is the incident angle [9].Therefore, a material's absorption characteristic can be completely understood if n and κ are known ( i θ is not a material property).Here, as well known [10], n and κ can be expressed in terms of primitive variables as ( ) where c is the speed of light in free space; ε , μ , and σ are permittivity, permeability, and electrical conductivity of the material; and λ is the laser wavelength in free space.If we look at Eqs. ( 4) and ( 5) carefully, we can learn that wavelength λ can be altered without changing n and κ if σλ remains unchanged.In other words, if we want to use a wavelength that is 100 times larger than the actual wavelength, we can decrease the electrical conductivity σ 100 times and still the angle-dependent absorption characteristic of the material remains exactly the same.Note that, in Eqs. ( 4) and (5) with n and κ values fixed, there are two equations in two unknowns, με and 2 c μσλ π , so these unknowns can be completely determined., i.e., Now, we can choose λ, ε, μ, and σ values by using Eq.(6).Note that in this study we let

Changing wavelength for dispersive FDTD algorithm
For many metals with high reflectance values, extinction coefficient κ is larger than refractive index n, and therefore, Eqs. ( 4) and ( 5) become inappropriate and the strategy presented in Section 2 cannot be employed.(Comparing Eqs. ( 4) and ( 5), everything is exactly the same except for the sign between the two terms in parentheses, so the latter cannot be larger than the former.)In this case, the standard FDTD method cannot be used to simulate metals, and the dispersive FDTD method needs to be used [6].In this section, we will present a scheme to use an enlarged wavelength, which can be used with the dispersive FDTD method.
The parameters for the dispersive FDTD are obtained from the Drude model [6] with the given n and κ values.Equations ( 4) and ( 5) can be re-written as where ε 1 and ε 2 are the real and imaginary parts of complex dielectric constant r ε and r μ is relative permeability.Also, the Drude model for the complex dielectric constant can be expressed as where ω p is plasma frequency, γ p is damping constant, and ω is the angular frequency of the laser beam, which is related to the beam wavelength as follows: Note that our objective is to increase λ (i.e., decrease ω) while maintaining the same n and κ values.From Eqs. ( 7)- (10), if 1 r μ = is assumed, we can increase λ as long as ε 1 and ε 2 are unchanged.In other words, from Eq. ( 9), when ω is changed, we can select proper values of ε ∞ , p ω , and p γ that will result in the same ε 1 and ε 2 values, which will in turn lead to the same n and κ values.

Using standard FDTD Scheme for n ≤ κ
In the previous two sections, we have discussed methods to increase wavelength for both standard and dispersive FDTD methods without changing the metal's absorption characteristics.The standard FDTD algorithm (Yee's original algorithm), however, can be employed only when n κ > (as is shown from Eqs. ( 4) and ( 5)), although it is simpler to implement than the dispersive FDTD algorithm, which can be used for virtually all n and κ values.In this section, we will present a scheme that enables the use of the standard FDTD scheme for n κ ≤ cases.In order to figure out the afore-mentioned issue when n κ ≤ , let's take a look at Eqs. ( 1) and ( 2).If the incident angle is fixed and the reflectance value is assigned, the two equations become quadratic equations in terms of n and κ.Here, for the sake of simplicity, the incident angle will be fixed at 0 o i θ = (i.e., at normal incidence), where p-and s-polarizations become identical.Figure 1 shows the reflectance value contour lines plotted using Eqs.( 1) and ( 2).In this figure, for a reflectance value, there exist infinitely many sets of n and κ , shown as a contour line.Furthermore, out of these infinite sets, we can notice that n κ > cases always exist.In other words, although n κ ≤ for the given material, we can always find a different set of n and κ with the same reflectance and n κ > .With a newly obtained set of n and κ with n κ > , now we can use the standard FDTD method instead of the dispersive FDTD method.
In fact, this result is true only when the incident angle is 0° in the whole computational domain.For simple problems like a flat plane vacuum-material interface, there is only one incident angle in the whole domain, and a single set of n and κ can be used to maintain the exact same reflectance according to Eqs. ( 1) and (2).In most problems, however, this is not the case and the incident angle cannot be assumed constant because the structure geometry could be very complicated: the incident angle can assume any values between 0° and 90°, and thus, different sets of n and κ values will be required at different locations and times, therefore, this scheme is virtually impossible to implement.
Then, naturally a question arises as to if there exists a certain representative incident angle * i θ that, when used in Eqs. ( 1) and ( 2), can resemble the original reflectance characteristic of the material very closely over the entire range of incident angle.To answer this question, let's consider an iron target irradiated by a 1.06 μm laser, where n = 3.81 and κ = 4.44 [9].In this case, apparently n κ ≤ , and we need to come up with a new set of n and κ, which reproduces the reflection characteristic of Fe at 1.06 μm over the entire range of incident angle.Figure 2 shows angle-dependent reflectance values of iron (Fe) under 1.06 μm irradiation.Here, the red curves represent the actual reflectance patterns of Fe constructed for s-and p-polarizations by using Eqs.( 1)-( 2) and n = 3.81 and κ = 4.44.In this case, the reflectance value at normal incidence is found to be 0.644.Now, let's choose * 0 i θ = °, and evaluate Eqs. ( 1) and (2) assuming 0.644 . Then, from the obtained quadratic equations we can take infinitely many sets of new n and κ, several of which are listed in Table 1.Using these new n and κ values, we can re-evaluate Eqs. ( 1) and ( 2) as a function of incident angle i θ , which are shown as blue lines in Fig. 2.
Here, we can notice two important things.First, for s-polarized lights, regardless of n and κ values, the calculated angle-dependent reflectance curves are very close to the actual reflectance curve of Fe at 1.06 μm irradiation (red curve) over the entire incident angle range from 0° to 90°.Secondly, for p-polarized lights, as the n value increases from the smallest to the largest, the reflectance curve approaches the actual reflectance curve from above and passes it eventually.Therefore, in between, there exist some n and κ values that very closely approximate the actual reflectance pattern of the material.Here, our strategy is to find a set of n and κ, which gives the most accurate approximation to the actual and also satisfies n κ > .For example, when n = 4.62 and κ = 4.51, the approximated reflectance curve is reasonably close to the actual especially if the incident angle is away from the minimum reflectance point.
, now we have two different quadratic equations for different polarizations.In order to estimate the error generated from this approximation, the following definition is used: ( ) Here, , ( )   s p R θ is the actual reflectance curve for p-or s-polarizations calculated by Eqs. ( 1) and ( 2), and * , ( ) R θ is the approximated reflectance curve with a new set of n and κ values.Figure 5 presents the relative errors for s-and p-polarizations as a function of the newly chosen n value.Here, we considered two * i θ values, 0° and 75.2°.As shown clearly, near the original value of n the error is small and it increases as it moves away from this point.The minimum errors occur at around n = 3.81 (original value), and this is understandable.Also, we can notice that when * 0 i θ = ° errors are smaller than when * 75.2 i θ = °.We can see that for s-polarized lights the error is very small for the almost entire range of n values while for ppolarization errors are much larger but still reasonably small if n is chosen near the original n value.For n = 4.62 and κ = 4.51, the relative errors for s-and p-polarizations are 0.04% and 3.4%, respectively.For p-polarization, the relative errors will be much smaller if the incident angle is not very large.Note that, although the scheme has been explained for a Fe target under 1.06 μm irradiation it can be used for other metal/wavelength combinations.

Results and discussion
For validation purposes, we have implemented the presented schemes in our in-house 3-D FDTD codes.We considered Sn targets under 1.06 μm irradiation (n = 4.7 and κ = 1.6 [9]) and Fe targets under 1.06 μm irradiation (n = 3.81 and κ = 4.44).For both materials, four types of simulations were performed: (a) standard FDTD simulation with the given wavelength, (b) dispersive FDTD simulation with the given wavelength, (c) standard FDTD simulation with an increased wavelength, and (d) dispersive FDTD simulation with an increased wavelength.In the case of the Sn target, n is already larger than κ, so a new set of n and κ does not have to be obtained for standard FDTD simulations ((a) and (c)).On the other hand, for the Fe target, a new set of n and κ needs to be used because n κ < .For all simulations, we calculated the laser beam absorptance as the incident angle changes from 0° to 70° at an interval of 10°.Note that 80° and 90° incident angles were not simulated because these require extremely large domain size to capture the entire reflection phenomena.Figure 6 shows the schematic diagram of the computational domain with dimensions when the original wavelength is used.When the wavelength is enlarged, we proportionally increased the domain size.In this study, the laser beam propagates in the negative z-direction, and the surface tilting angle θ is increased from 0° to 70°.For all cases, p-, s-, and circularly polarized Gaussian beams were considered.Table 2 and Table 4 present the parameters used for the standard FDTD simulation of Sn and Fe, respectively, and in Table 3 and Table 5 are presented the parameters for the dispersive FDTD simulations of Sn and Fe, respectively.All these parameters were obtained by using the methods presented in Sections 2 to 4. For the Sn target, the wavelength has been increased to 31.8 μm (30 times larger), and we have used a wavelength of 21.2 μm (20 times larger) for the Fe target.These wavelengths were selected arbitrarily, and of course, much larger wavelengths can be used as long as the electrical conductivity is large enough to absorb the electromagnetic waves near the material interface.For all simulations, uniform grids with a grid spacing of / (10 )  n λ was used, and a 12 core high performance computer was used.Running times varied from 4 to 8 hours depending on the problem size.Figure 7 shows the simulation results showing electric fields for different tilting angles.In this simulation, the dispersive FDTD algorithm was used for a Fe target with an enlarged wavelength and the beam was s-polarized.In this figure, we can clearly see how the beam reflection changes as the tilting angle increases.Note that the tilting angle is equal to the incident angle, and when it is larger than 50°, the z dimension has to be increased to capture the whole reflection phenomena.For an incident angle of 80°, as mentioned earlier, the required domain becomes extremely large, so that case was omitted.However, the authors believe that the simulation results for larger angles will be qualitatively similar.
In Fig. 8, the calculated reflectance versus incident angle for all simulations is presented.In this study, the reflectance was calculated by using the energy flux difference between the incident waves and reflected waves in terms of the Poynting vector.As shown clearly in the figures, for both Sn and Fe, for all polarizations, and for all simulation methods, the obtained results are in good agreement with the analytical solutions.Considering that uniform grids were used and the surface is not aligned with any of the grid lines, we believe that the simulation results are reasonably accurate.In order to validate the method with a more complex geometry, we considered a three dimensional problem, where a 1.06 μm Gaussian beam irradiates on the upper side of a cylinder made of Fe. Figure 9 shows the schematic diagram of the computational domain with dimensions when the original wavelength is used.The same problem was also solved with a 20 times larger wavelength (21.2 μm), and in this case all the dimensions were increased proportionally.In this study, we have conducted four types of simulations for validation purposes: (a) standard FDTD simulation (new n and κ) with λ = 1.06 μm, (b) dispersive FDTD simulation with λ = 1.06 μm, (c) standard FDTD simulation (new n and κ) with an increased wavelength of 21.2 μm, and (d) dispersive FDTD simulation with an increased wavelength of 21.2 μm.For each case, three different polarizations were considered: p-, s-, and circular polarizations.Apparently, if the electric field of the beam is aligned in the y- direction, even though the cylinder surface is curved, the laser beam is 100% s-polarized.Also, if the electric field is in the x-direction, the laser beam is entirely p-polarized.All simulation parameters are listed in Table 4 and Table 5.  Figure 10 shows the simulation results for Case (d).It turned out that all other results are undistinguishably similar to this result, so in this study only one result is shown.In the first and second rows, the electric fields on the x-y and y-z planes are presented, respectively, where both planes pass through the center of the laser beam.For all three polarization results, the electric fields are astonishingly similar, and the laser beams look smaller when viewed through the x-z plane because of more complex interference of incident and reflected waves.To validate the results, we have analytically calculated the reflectance values of the Fe cylinder for s-and p-polarizations as follows: where ℜ and 0 r are cylinder and beam radii, respectively, and 0 E is the electric field value at the center of the Gaussian beam.Note that the reflectance for circular polarization is mathematically the average of * s R and * p R .Table 6 presents the simulated reflectance values and the corresponding errors with respect to the analytical solutions obtained from Eq. (12).As we can see from the results, all four cases agree with the analytic solutions well, and the maximum relative error is less than 4.6%.Especially, it is clearly shown from the table that the result obtained with an increased wavelength is virtually the same as the one from the corresponding original-wavelength simulation.This means that the relative errors shown in the table (up to 4.6%) are the errors of the FDTD method, and not the errors caused by the additional procedures proposed in this study.Moreover, we can notice that the dispersive FDTD simulations are found to be slightly more accurate than the standard FDTD ones, which is because newly selected n and κ values were used for the standard FDTD simulations of Fe.For dispersive and standard FDTD simulations, the maximum errors are 3.27% and 4.59%, respectively, both of which occurred for p-polarization.Also, from the results, we can see that the reflectance for circular polarization is exactly the average of the p-and spolarizations.
One last comment is that, if an enlarged wavelength is used, the actual wave characteristic such as diffraction and interference are also changed although the beam absorption characteristic is preserved.However, in many problems, such as in laser processing problems, the most important things are the laser beam absorptivity and the beam absorption pattern.Besides, in laser manufacturing problems, the focused laser beam diameter is ~500 μm (which is much larger than the wavelength), so the beam divergence is generally very small compared to the problem size.

Conclusions
In this article, we have presented and validated an FDTD based method for simulating largescale laser beam absorption problems by enlarging laser wavelength while maintaining angledependent absorption characteristics.A method to use the standard Yee algorithm for a material with n κ < has bee also presented.Using these methods, we believe that various problems where laser beam absorption is critical can be effectively solved.

Table 1 .
Several selected n and κ values that lead to the same reflectance at normal incidence (θ i * = 0°) as the actual reflectance value of iron under 1.06 μm irradiation (n = 3.81, κ = 4.44).Bold faced cases are shown in Fig. 2 as blue lines.

Figure 3 (
b) shows the contours of a reflectance value of 0.644 for s-and ppolarizations.The original n and κ values of Fe are shown at the intersection of two curves as a green circle.In finding a new set of n and κ values that satisfies n κ > , even though it is always possible, we notice one problem.Because now the two polarization cases have different quadratic curves, for different polarizations a different set of n and κ values need to be used.If a light polarization is well defined and fixed in a given problem, it is not a problem at all.However, in most problems, light polarization is arbitrary and/or changes from one place to another, so that having to select different n and κ values for different polarizations is impractical.Furthermore, as shown in Fig.4, the calculated reflectance curves, especially for p-polarization, is much worse except when * 75.2 i θ = °.In this study, we have also tested other * i θ values (30°, 45°, and 82°), but * 0 i θ = ° were found to be the most desirable because of the same reasons explained above.

Fig. 3 .Fig. 4 .
Fig. 3. Lines of constant reflectance (R = 0.644) for s-and p-polarized lights shown as blue dashed lines and red solid lines, respectively.

Fig. 5 .
Fig. 5. Overall relative errors generated when a different n and κ values are used to approximate the material's reflectance patterns.Errors are calculated for Fe under 1.06 μm irradiation.

Fig. 7 .
Fig. 7. Electric fields in the computational domain showing the reflection patterns at different tilting angles.Simulations were performed for Fe with an enlarged wavelength using the dispersive FDTD code assuming the laser beam is s-polarized.The width of the computational domain is 200 μm.

Fig. 10 .
Fig. 10.Simulation results showing a Gaussian beam irradiating on a cylinder made of Fe.Here, the dispersive FDTD with an enlarged wavelength of 21.2 μm was used instead of 1.06 μm.The width and height of the figures are both 240 μm.
was totally arbitrary, and our goal is to obtain the best approximation.Therefore, we need to know what * i θ gives the best result.From the above result, because the error is the largest near the minimum reflectance point, we will try to use Brewster's angle B θ , which is 75.2° in this case.If Eqs.(1) and (2) are evaluated using i θ = °