Printer Friendly

Hybrid Simulation of Maxwell-Schrodinger Equations for Multi-Physics Problems Characterized by Anharmonic Electrostatic Potential.


Hybrid simulation is an essential approach to study multi-physics phenomena ruled by more than one governing equations. One of typical multi-physics examples is the light-matter interaction in which an incident laser field excites electrons in the matter and then the excited electrons in turn modify the electromagnetic field as a polarization current source. In the last decade this cooperative light-matter interaction has been actively studied with growing interests particularly in designing innovative photonic devices such as plasmonic antennas, quantum dot laser, and so on [1-3]. Currently, the multi-physics phenomena in such devices have been studied by two distinct hybrid schemes, namely, based on the Maxwell-Schrodinger and Maxwell-Newton theories [4-9]. The Maxwell-Schrodinger hybrid scheme [4-8] is computationally demanding but physically precise, where the laser field is governed by Maxwell's equations and the electrons by Schrodinger's equation. This scheme has been successfully used recently in pioneering numerical simulations such as for a carbon nanotube transistor [4], [H.sup.+.sub.2] gas interacting with ultrashort laser pulses [7,8], and so on. The other approach, namely, the well-known Maxwell-Newton hybrid scheme [9], where the laser field and electrons are described by solving Maxwell's equations and Newton's equation of motion, respectively, requires much less computational resources, and thus has been very often employed to perform multi-physics simulations [1-9] without, however, paying much attention to its physical reliability. Very recently, we have examined a reliability of the conventional Maxwell-Newton hybrid scheme for an electron confined in one-dimensional potential wells. A comparison of computational results with those obtained by the corresponding Maxwell-Schrodinger scheme [6] has revealed that both results agree very well for a harmonic single-well potential while disagree qualitatively for a triple-well potential where a harmonic potential is artificially supplemented by two small humps allowing bifurcation of the electron wave packet. Although this study clearly demonstrates necessity of use of the Maxwell-Schrodinger scheme over the conventional Maxwell-Newton one particularly when quantum-mechanical tunnelling takes place, an effect of anharmoicity in a single-well confining potential on computational results has not been fully explored.

In this paper we have studied a system of a nano-scale thin film interacting with pulsed laser fields whose electrostatic confining potential for electron is modelled by a single-well potential. Since tunnelling effects can be safely neglected for a single-well potential, we can make unambiguous manifestation of an effect of anharmonicity in the confining potential on differences between computational results obtained from the two distinct hybrid schemes. This allows us to clarify further the extent of applicability of the conventional Maxwell-Newton scheme from physical viewpoints.


Figure 1 illustrates our theoretical model used in the present study. The thin film is uniform in the y-z plane and its optical properties are assumed to be calculated from the responses of one representative electron among a larger number of electrons comparable to the order of Avogadro's number. The incident electromagnetic fields consisting of only [E.sub.y] and [H.sub.y] components are given by a plane wave, which excite all electrons in the film to the polarization direction y. This one-dimensional model enables us to solve both of Maxwell-Schrodinger and Maxwell-Newton equations very accurately, allowing us a detailed comparison of their computational results.

2.1. Maxwell-Schrodinger Hybrid Scheme

The computational procedure to solve the Maxwell-Schrodinger scheme is described in Figure 2(a). Maxwell's equations for dielectric objects are given by

[nabla] x E = [[mu].sub.0] [partial derivative]H/[partial derivative]t, (1)

[nabla] x H = [[epsilon].sub.0] [partial derivative]H/[partial derivative]t + J, (2)

where J represents the polarization current density which is defined by the time derivative of the polarization vector P.

Since the electromagnetic fields have only [E.sub.y] and [E.sub.z] components in the present study, they can be updated by the following recursion relations based on the Maxwell FDTD method [5, 6,10]:

[mathematical expression not reproducible], (3)

[mathematical expression not reproducible], (4)

where n, i, iF, and [delta] represents the time step, the space grid along the x axis, the cell position of the thin film, and Kronecker delta function, respectively. The edges of computational domain in the Maxwell FDTD simulation are supplemented by the Mur absorbing boundary condition [10].

The Schrodinger equation for an electron subjected to a laser field is given by

[mathematical expression not reproducible], (5)

where the so-called length gauge has been adopted to describe the interaction between the electron and the electromagnetic field [11]. The following recursion relations based on the Schrodinger FDTD method [4-8,12] can be obtained by separating the real and imaginary parts of the Schrodinger equation:

[mathematical expression not reproducible], (6)

[mathematical expression not reproducible], (7)

where [[psi].sub.mag] and [[psi].sub.real] are the imaginary and real parts of the wave function [psi] discretized on the space grids {j} placed along the y axis. The operator a in these equations performs the following sixth-order accurate difference to simulate the second-order derivative [[partial derivative].sup.2]/[partial derivative] [y.sup.2] for an arbitrary function F:

[mathematical expression not reproducible]. (8)

We employ the following Dirichlet boundary for the Schrodinger FDTD simulation:

[mathematical expression not reproducible], (9)

where [j.sub.max] denotes the number of grids for the y axis. This condition is well-known to induce spurious oscillations when the wave function impinges on the boundary. Therefore, we utilize sufficient wide analysis domain so as to avoid these numerical artifacts.

The polarization current density J in the Maxwell-Schrodinger scheme is defined by the following expression with the electron density N:

J = qN [[integral].sup.[infinity].sub.-[infinity]] [[psi].sup.*] [nabla] [psi] dv. (10)

Equation (10) describes the average behaviour of the current density due to the motion of all electrons expressed by the wave function of a representative electron. The y component of (10) can be evaluated by

[mathematical expression not reproducible]. (11)

where the operator [beta] performs the following sixth-order accurate difference to simulate the first-order derivative [partial derivative]/[partial derivative]y for an arbitrary function F.

[mathematical expression not reproducible]. (12)

The Maxwell-Schrodinger scheme is realized by using Equations (3), (4), (6), (7), and (11) recursively as illustrated in Figure 2(a).

2.2. Maxwell-Newton Hybrid Scheme

The computational procedure adopted in the Maxwell-Newton scheme is shown in Figure 2(b). The part for solving Maxwell's equations is the same as in the Maxwell-Schrodinger schemes based on (3) and (4). The following Newton equation is employed to describe the motion of an electron confined by the electrostatic potential V and subjected to an external electromagnetic field:

m [d.sup.2]/[dt.sup.2] = q[E.sub.y] + Fv, (13)

Fv = -[partial derivative]V/[partial derivative]y. (14)

where we assume that the electron feels no frictional force. The polarization vector P and polarization current density J in the Maxwell-Newton scheme [9] are, respectively, defined by

P = qNr, (15)

J = [partial derivative]P/[partial derivative]t. (16)

One can derive the following recursion relations for simulating these polarization and current density in the FDTD framework as

[mathematical expression not reproducible], (17)

[mathematical expression not reproducible]. (18)

In the Maxwell-Newton scheme Equations (3), (4), (17), and (18) are solved recursively as displayed schematically in Figure 2(b).


The incident laser fields are generated from the following electric and magnetic current sources [J.sup.(i).sub.e] and [J.sup.(i).sub.m] with the unit function u(t)

[mathematical expression not reproducible], (19)

[mathematical expression not reproducible], (20)

where we have set [J.sub.0], [DELTA]x, [sigma]t, and [t.sub.0] as 1000 MA/m, 0.125 nm, 1.25 fs, and 20[[sigma].sub.t] [[sigma].sub.t], respectively. The time step [DELTA]t is chosen to be smaller by a factor of 0.9 than [DELTA][t.sub.CFL], i.e., the maximum value allowed in the CFL condition [10], so as to guarantee numerical stability.

We compare the results simulated by the Maxwell-Schrodinger and Maxwell-Newton schemes for the following three electrostatic potentials [V.sub.S], [V.sub.D], and [V.sub.AH]:

[V.sub.S] = m[w.sub.s][y.sup.2], (21)

[V.sub.D] = [V.sub.s] + [V.sub.0] exp {-0.5 ([y.sub.D] + y/[[sigma].sub.D]).sup.2]}, (22)

[V.sub.AH] = [V.sub.1] [(y/[y.sub.AH]).sup.4], (23)

where the parameters characterizing the potentials, [w.sub.S], [V.sub.0], [[sigma].sub.D], [y.sub.D], [V.sub.1], and - [y.sub.AH], are given as 50Trad/s, 0.5 eV, 0.625 nm, 10[[sigma].sub.D] nm, 4.5 eV, and 25 nm, respectively. The potential energy curves for these three potentials are plotted in Figure 3. As displayed in this figure Vs is a single-well and harmonic potential, while [V.sub.D] is almost identical to this [V.sub.S] potential but is supplemented by a small 'humps' located at around - = -6.25 nm. This hump allows the quantum electron to bifurcate every time when it impinges on the hump owing to tunnelling while does not for the classical electron, as we explored in our previous study [6]. The third potential [V.sub.AH] is single-well but anharmonic, which allows us to investigate a different quantum mechanical effect other than that caused by tunnelling. We have chosen the ground state of each of these electrostatic potentials as the initial wave packet in all quantum simulations.

The time responses of the polarization current density J in the thin film obtained by the two hybrid simulations are represented in Figure 4, where the blue solid and red broken lines represent, respectively, the numerical results obtained by the Maxwell-Schrodinger and Maxwell-Newton schemes. Figure 4(a) representing the results for the single and harmonic well Vs shows that both results agree excellently, indicating that the classical theory of the Maxwell-Newton scheme can be safely used for this case. On the other hand, Figure 4(b), displaying the results for the double-well potential VD, shows that the polarization current densities obtained from these two schemes deviates from each other more and more strongly after the first 30 fs. This indicates that the Maxwell-Newton scheme is unreliable for this double-well potential, as was demonstrated in our previous study for a triple-well potential [6]. The results displayed in Figure 4(c) for the anhamonic potential [V.sub.AH] shows a trend somewhat between (a) and (b): the polarization current density of the Maxwell-Newton scheme roughly follows that of the Maxwell-Schrodinger scheme, but there can be observed a quantitative difference between them. This indicates that the Maxwell-Newton scheme could become unreliable for quantitative calculation even when the confining potential is single well.

In order to rationalize the observed trends, we have investigated the dynamics of electron in the thin film, namely, the spatiotemporal propagation of the electron wave packets and the corresponding classical trajectories obtained, respectively, by the Maxwell-Schrodinger and Maxwell-Newton schemes. The results for the three electrostatic potentials [V.sub.S], [V.sub.D], and [V.sub.AH] are displayed in Figures 5(a), 5(b), and 5(c), respectively. In the figures the thick oscillatory curve in color whose scale is displayed on the right end of each figure indicates the time-evolution of the probability density of the electron wave packet [[absolute value of [psi]].sup.2] and the triangles plotted in the same figure denote the classical trajectory. On the left-hand side of each figure the potential energy curve of the corresponding electrostatic potential is also plotted. The vertical axes for both sides of the figure commonly indicate the y axis. As shown in Figure 5(a) representing the results for the single and harmonic well, the electron wave packet is localized at each time step keeping a Gaussian shape similar to the ground state and closely follows the corresponding classical trajectory. This excellent agreement between the quantum and classical electron dynamics results in an almost identical behaviour of the current densities obtained by these two schemes as displayed in Figure 4(a). On the other hand, Figure 5(b), representing the results for the double-well [V.sub.D], shows that the electron wave packet gets fragmented into several pieces due to tunnelling. Furthermore interference among these fragments makes the wave packet complicated even further. Since the classical dynamics could not support such fragmentation and interference, the current density obtained by the Maxwell-Newton scheme deviates largely from that obtained by the Maxwell-Schrodinger scheme as observed in Figure 4(b). Figure 5(c), representing the results for the anharmonic potential [V.sub.AH], shows that the electron wave packet follows the corresponding classical trajectory in the beginning before t~200fs. For the later time t, however, the electron wave packet starts to spread gradually and a nodal structure in the probability density appears. This nodal structure reflects the fact that the electron wave packet is no more a single Gaussian distribution but is fragmented into a few components. In case for purely harmonic electrostatic potentials an initial Gaussian wave packet remains to be a Gaussian through time propagation. Therefore, the observed fragmentation is caused by anharmonicity in the electrostatic potential, which induces dephasing of the electron wave packet. Since classical mechanics cannot account for such dephasing effects, the classical trajectory deviates from the center of the electron wave packet, which causes a difference in the polarization current density between the two schemes.

Next, we have examined a dependence of the computational results on the strength of the applied laser field. Since the applied laser field we have studied so far is rather strong, we have employed a weaker laser field here by decreasing the amplitude of the current sources 10 times smaller than that used for the simulations in Figures 4 and 5 as [J.sub.0] = 100MA/m. Figures 6(a) and 6(b) display the resultant time responses of the polarization current densities J for the single- and double-well potentials [V.sub.S] and [V.sub.D]. The blue solid and red broken lines represent the results obtained by the Maxwell-Schrodinger and Maxwell-Newton schemes, respectively, as for Figure 4. Unlike the results in Figure 4 the polarization current densities for not only the single-well potential Vs but also the double-well potential [V.sub.D] obtained by the two schemes agree very well as displayed in Figures 5(a) and 5(b). On the other hand, Figure 6(c) representing the result for the anharmonic single-well potential [V.sub.AH] shows that the computational results by the two hybrid simulations still differs quantitatively from each other.

As has been done in Figure 5 for the strong electromagnetic field, the time evolution of the electron wave packet and the corresponding classical trajectory for this weak laser field are displayed in Figure 7, where (a), (b), and (c) denote the numerical results for [V.sub.S], [V.sub.D], and [V.sub.AH], respectively. It is noted that Figures 7(a) and 7(b), representing the results for [V.sub.S], [V.sub.D], are almost identical to each other and that no fragmentation of the electron wave packet is observed for the double-well case unlike the corresponding result in Figure 5(b). This can be rationalized by the small strength of the laser field as follows: since the electric field of the laser pulse is small, it could not give enough energy to the electron to reach the hump of the potential for [V.sub.D]. Therefore, since the potential energy curves of Vs and [V.sub.D] below this hump are exactly the same harmonic potential, their electron dynamics should naturally be identical to each other for this weak strength of the laser field. In the case of the anharmonic potential [V.sub.AH] illustrated in Figure 7(c), however, the electron wave packet spreads as the time proceeds and it undergoes bifurcation after t = 300 fs. Therefore, this dephasing effect existing only in the quantum simulation causes a difference between the results obtained by the two schemes even when the laser field is sufficiently weak.

The present investigations show that the conventional Maxwell-Newton scheme can be applied not only to macroscopic problems as have been studied in most cases but also to microscopic problems of a nano-scale order on condition that the electrostatic confining potential for electron is purely harmonic. However, when the electrostatic potential deviates from a harmonic one even slightly, the Maxwell-Newton scheme would give unreliable results owing to quantum-mechanical tunnelling and/or anharmonicity effects. Therefore, such problems should be solved by the Maxwell-Schrodinger hybrid scheme.


In this paper, we have focused on the anharmonisity of the electrostatic potential, and investigated the interaction between laser fields and a nano-scale thin film modelled by an electron confined in an electrostatic potential. The two distinct hybrid simulations, the Maxwell-Schrodinger and the conventional Maxwell-Newton schemes, have been compared for an anharmonic single-well potential problem, where quantum-mechanical tunnelling does not take place. Furthermore, harmonic singlewell and harmonic double-well potential problems have been also investigated to make a comparison with the results for the anharmonic single-well potential. The computational results show that the two multi-physics simulations provide almost identical results for the harmonic confining potential, indicating a validity of use of the conventional Maxwell-Newton scheme for this case. In the case of the double-well potential, however, the results by the Maxwell-Newton approach differ significantly from those by the Maxwell-Schrodinger approach when tunnelling plays an important role. Finally, for the case of the anharmonic potential, the result of the Maxwell-Newton simulation deviates from that of the Maxwell-Schrodinger simulation quantitatively owing to an effect of dephasing of the electron wave packet by the anharmonicity in the electrostatic potential. These results have clearly demonstrated that the Maxwell-Schrodinger scheme is indispensable to multi-physics simulation particularly when the tunnelling, interference and anharmonicity effects play an essential role.


The authors would like to thank Prof. K. Nakagawa, Prof. Y. Ashizawa (Nihon University), Prof. M. Tanaka (Gifu University), and Prof. W. C. Chew (University of Illinois) for their useful comments and suggestions. This work was partly supported by Grant-in-Aid for Scientific Research (C) (26420321) and MEXT-Supported Program for the Strategic Research Foundation at Private Universities, 2013-2017.


[1.] Ota, T., Y. Ashizawa, K. Nakagawa, S. Ohnuki, H. Iwamatsu, A. Tsukamoto, and A. Itoh, "Dependence of circularly polarized light excited by plasmon aperture on relative position to magnetic particles for all-optical magnetic recording," Journal of the Magnetics Society of Japan, Vol. 36, Nos. 1-2, 66-69, Feb. 2012.

[2.] Sugawara, M., N. Hatori, M. Ishida, H. Ebe, Y. Arakawa, T. Akiyama, K. Otsubo, T. Yamamoto, and Y. Nakata, "Recent progress in self-assembled quantum-dot optical devices for optical telecommunication: Temperature-insensitive 10Gbs_1 directly modulated lasers and 40Gbs_1 signal-regenerative amplifiers," Journal of Physics D: Applied Physics, Vol. 38, No. 13, 2126-2134, Jun. 2005.

[3.] Aoki, K., D. Guimard, M. Nishioka, M. Nomura, S. Iwamoto, and Y. Arakawa, "Coupling of quantum-dot light emission with a three-dimensional photoniccrystal nanocavity," Nature Photonics, Vol. 2, No. 13, 688-692, Oct. 2008.

[4.] Pierantoni, L., D. Mencarelli, and T. Rozzi, "A new 3-D transmission line matrix scheme for the combined Schrodinger-Maxwell problem in the electronic/electromagnetic characterization of nanodevices," IEEE Transactions on Microwave Theory and Techniques, Vol. 56, No. 3, 654-662, Mar. 2008.

[5.] Ohnuki, S., T. Takeuchi, T. Sako, Y. Ashizawa, K. Nakagawa, and M. Tanaka, "Coupled analysis of Maxwell-Schroodinger equations by using the length gauge: Harmonic model of a nanoplate subjected to a 2D electromagnetic field," International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, Vol. 26, 533-544, Feb. 2013.

[6.] Takeuchi, T., S. Ohnuki, and T. Sako, "Comparison between Maxwell-Schrodinger and Maxwell-Newton hybrid simulations for multi-well electrostatic potential," IEEE Journal of Quantum Electronics, Vol. 50, No. 5, 334-339, May 2014.

[7.] Lorin, E., S. Chelkowski, and A. D. Bandrauk, "A numerical Maxwell-Schrodinger model for intense laser-matter interaction and propagation," Computer Physics Communications, Vol. 177, No. 12, 908-932, Jul. 2007.

[8.] Lorin, E., S. Chelkowski, and A. D. Bandrauk, "Attosecond pulse generation from aligned molecules-dynamics and propagation in H+," New Journal of Physics, Vol. 10, No. 2, Feb. 2008.

[9.] Yamaguchi, T. and T. Hinata, "Optical near-field analysis of spherical metals: Application of the FDTD method combined with the ADE method," Optics Express, Vol. 15, No. 18, 11481-11491, Sep. 2007.

[10.] Taflove, A. and S. C. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, 3rd edition, Artech House, Boston, London, 2005.

[11.] Cohen-Tannoudji, C., J. Dupont-Roc, and C. Crynberg, Atom-photon Interactions: Basic Processes and Applications, Wiley-VCH, Weinheim, 2004.

[12.] Soriano, A., E. A. Navarro, J. A. Porti, and V. Such, "Analysis of the finite difference time domain technique to solve the Schrodinger equation for quantum devices," Journal of Applied Physics, Vol. 95, No. 12, 8011-8018, Jun. 2004.

Takashi Takeuchi, Shinichiro Ohnuki *, and Tokuei Sako

(Invited Paper)

Received 30 June 2014, Accepted 19 July 2014, Scheduled 21 July 2014

Invited paper for the Commemorative Collection on the 150-Year Anniversary of Maxwell's Equations.

* Corresponding author: Shinichiro Ohnuki (

The authors are with the College of Science and Technology, Nihon University, 1-8-14 Surugadai, Kanda, Chiyoda-ku, Tokyo 101-8308, Japan.

Caption: Figure 1. Geometry and coordinate systems. A thin film and incident current sources, illustrated by a grey box and blue arrows, respectively, are uniform in the y-z plane. All electrons in the film are confined in the electrostatic potential V, and can move along the y axis which is parallel to the direction of the electric field.

Caption: Figure 2. A schematic illustration of the computational schemes for the two hybrid simulations: (a) Maxwell-Newton and (b) Maxwell-Schrodinger.

Caption: Figure 3. Spatial profile of the studied electrostatic potentials: the blue line with circles represents the harmonic single-well potential Vs while the red and green lines denote the almost harmonic double-well potential [V.sub.D] and anharmonic single-well potential [V.sub.AH], respectively.

Caption: Figure 4. Comparison of the time response of the polarization current density for the electrostatic potentials [V.sub.S], [V.sub.D], and [V.sub.AH] (See Figure 3). (a), (b), and (c) Correspond, respectively, to the case for the electrostatic potential [V.sub.S], [V.sub.D], and [V.sub.AH]. The blue solid and red broken lines represent the results obtained by the Maxwell-Schrodinger and Maxwell-Newton schemes, respectively.

Caption: Figure 5. Time evolution of the electron wave packet and the corresponding classical trajectory. (a), (b), and (c) Correspond, respectively, to the results for the electrostatic potential [V.sub.S], [V.sub.D], and [V.sub.AH]. The spatial profile of the potential is displayed on the left-hand side of each figure. The thick curve in color scale represents the probability density of the electron.

Caption: Figure 6. Comparison of the time response of the polarization current density for different electrostatic potentials. (a), (b), and (c) Correspond, respectively, to the cases of [V.sub.S], [V.sub.D], and [V.sub.AH]. The thin film is subjected to the weak electromagnetic fields excited by current sources with the amplitude Jo = 100MA/m. See the caption to Figure 4 for other remarks.

Caption: Figure 7. Time evolution of the electron wave packet and the corresponding classical trajectory for the weak electromagnetic fields excited by current sources with the amplitude [J.sub.0] = 100MA/m. (a), (b), and (c) Correspond, respectively, to the cases of the electrostatic potential [V.sub.S], [V.sub.D], and [V.sub.AH]. See the caption to Figure 5 for other remarks.
COPYRIGHT 2014 Electromagnetics Academy
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2014 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Takeuchi, Takashi; Ohnuki, Shinichiro; Sako, Tokuei
Publication:Progress In Electromagnetics Research
Date:May 1, 2014
Previous Article:A Printed Vivaldi Antenna with Improved Radiation Patterns by Using Two Pairs of Eye-Shaped Slots for UWB Applications.
Next Article:Differential Forms and Electromagnetic Field Theory.

Terms of use | Privacy policy | Copyright © 2022 Farlex, Inc. | Feedback | For webmasters |