text:
# Scaling relations in early-type galaxies belonging to groups

## 1 Introduction

Galaxy scaling relations, such as the Fundamental Plane, and its various projections, and the colour-magnitude relation have been used to enhance our understanding of galaxy structure and evolution. Early-type galaxies, in particular, form a relatively homogeneous population and hence rather well-defined scaling relations. Most large samples of early-type galaxies come from the rich environments of galaxies clusters. Studies of these high density regions benefit from the large numbers of galaxies in a small angular region of the sky (hence making observations more efficient) and from the morphology-density relation, which tells us that such environments are dominated by early-type galaxies. Thus our current knowledge gained from galaxy scaling relations applies mostly to clusters. The field and group environments are less well studied in this respect, and provide a means to study environmentally-dependent processes. For example, galaxy mergers (which can be about hundred times more efficient in today’s groups than clusters, cf. Mamon 1986) may result in merger-remnant Ellipticals that deviate strongly from the scaling relations (e.g. Forbes, Ponman, & Brown, 1998). On the other hand many of the well studied Elliptical galaxies are brightest group galaxies.

Recent X-ray observations of groups and their constituent galaxies have provided a new insight into the study of the evolution of galaxies in small ensembles and in intermediate-density environments (Helsdon & Ponman, 2000; Mulchaey, 2000). Miles et al. (2003) find that the shape of the luminosity function of galaxies in groups with low X-ray luminosity is significantly bimodal, in comparison with that of X-ray bright groups or rich clusters, showing a strong environmental effect on the morphological mix in a system, and on the evolution of galaxies in general. It is likely that the dip in the luminosity function at intermediate magnitudes arises from increased merging in low velocity dispersion systems.

Of particular interest in this context are early-type galaxies, which are more abundant in clusters and groups than in the field. Several well-known scaling relations (such as the Faber-Jackson relation or the Fundamental Plane of Elliptical galaxies) utilise galaxy velocity dispersions and hence require spectroscopic data. However, scaling relations like the Kormendy relation (Kormendy, 1977) and the Photometric Plane of galaxies (Khosroshahi et al., 2000a) can be constructed from photometric observations alone.

There have been many studies aimed at understanding the differences in the structure of Ellipticals (E) and dwarf Ellipticals (dE), with mixed results. While some (e.g. Terlevich, Caldwell, & Bower, 2001) argue in favour of the similarities between the two populations, others (e.g. Kormendy, 1985) find evidence otherwise (see Graham & Guzmán 2003 for more details).

Traditionally, the surface brightness distribution of Ellipticals and bulges of disk galaxies have been modelled by the de Vaucouleurs profile. However, during the last decade, better photometry and advanced analysis have shown that their surface brightness distribution is, in general, better represented by the Sersic profile ($`r^{1/n}`$, Sersic 1968), of which the de Vaucouleurs profile ($`r^{1/4}`$) is a special case (Caon et al., 1993; Khosroshahi et al., 2000a; Graham, 2001). The Sersic model is also widely used to describe the brightness distribution of bulges of disk galaxies (de Jong, 1996; Andredakis et al., 1995; Khosroshahi et al., 2000b). The disk component of galaxies is usually characterised by an exponential distribution.

In this paper we examine the surface brightness distribution of Ellipticals belonging to galaxy groups, where most galaxies in the Universe reside. The relatively small number of galaxies in individual groups, compared to clusters, requires any analysis to be carried out on a large sample. This can be done in several ways. Here, the galaxies are classified based on the X-ray luminosity ($`L_X`$) of their parent groups, which, being correlated with the velocity dispersion of the group (e.g. Helsdon & Ponman, 2000), can be an index of the dynamical state and the mass of the group (Miles et al., 2003; Osmond & Ponman, 2003). The principal properties of the groups, and a detailed description of the sample can be found in section 2. The analysis techniques and morphology of galaxies is discussed in section 3. Correlations among various parameters and the scaling relations is studied in section 4. Section 5 contains a discussion and a summary of results.

## 2 Observations and Analysis

### 2.1 Observations

Our sample is drawn from the Group Evolution Multi-wavelength Study (GEMS, Osmond & Ponman, 2003) of sixty groups, compiled by cross-correlating a list of over 4000 known groups from various studies with the list of archival ROSAT PSPC X-ray observations with integration $`>10`$ ks. This includes several Hickson groups, as well as loose groups from the CfA survey. A large fraction of these were detected in the X-ray, and for the others we have upper limits for their X-ray luminosity.

Of the GEMS sample, 16 groups were observed at the 2.5m Isaac Newton telescope at the Roque de Los Muchachos Observatory, La Palma, between 2000 February 4–10. This is a random and representative subsample of the original sample of 60 groups, where all groups accessible for the allocated observing run were observed. The detector used was the Wide Field Camera (WFC), which is an array of four thinned EEV CCDs situated at the f/3.29 prime focus of the INT, each with an area of $`2048\times 4096`$ pixels, each pixel being 0.33 arcsec across. Each CCD thus can image an area of $`22.5\times 11.3`$ arcmin of sky, together covering 1017 sq.arcmin.

Photometry was performed with broadband BRI filters, of which we analyse only the $`R`$-band images here. Our analysis is limited to galaxies brighter than $`M_R=13`$. The median seeing achieved was about 1.1 arcsec in $`R`$. Further details can be found in (Miles et al., 2003).

### 2.2 Group Membership

The identification of the group galaxies is based on a colour selection. Source extraction was performed using the Sextractor software Bertin & Arnouts (1996), which uses a routine based on neural network generated weights for star-galaxy separation. Detections were checked visually and objects with the stellaricity parameter greater then 0.9 were deemed to be definitely stellar and therefore discarded. To obtain colours, a fixed aperture, set to be slightly greater than the seeing, was used to evaluate fluxes with different filters. All objects with FWHM less than the mean seeing were not processed further.

Galaxies were selected as being likely group members on the basis of their $`BR`$ colour. A conservative cut of $`BR=1.7`$ was found to remove the majority of background galaxies for all groups in our sample (See Miles et al. 2003 for further details). Those galaxies found in the NASA/IPAC Extragalactic Database (NED) within the virial radius $`R_{\mathrm{vir}}`$ (see Osmond & Ponman, 2003).

### 2.3 Analysis

The core of our data analysis is the two-dimensional decomposition of galaxies, by fitting a Sersic profile to the bulge and an exponential profile to the disk component, for obtaining the structural photometric parameters, and morphological classification. This analysis was carried out on all galaxies identified as group members using the selection method described in the previous section.

The bulge-disk decomposition procedure used here involves fitting two-dimensional image models to the observed galaxy images, assumed to consist of a bulge and a disk component. The model image is convolved with a Gaussian PSF, with the FWHM measured from the galaxy frame. Details of the accuracy and reliability of the decomposition procedure, and the associated galaxy simulation code, can be found in Wadadekar et al. (1999).

The Sersic profile for the bulge is given by

$`I_{\mathrm{bulge}}(x,y)`$ $`=`$ $`I_b(0)\mathrm{exp}\left[2.303b_n(r_{\mathrm{bulge}}/r_e)^{1/n}\right];`$ (1)
$`r_{\mathrm{bulge}}`$ $`=`$ $`\sqrt{x^2+y^2/(1ϵ_b)^2}.`$
Here $`x`$ and $`y`$ are the distances from the centre of the galaxy along the major and minor axis respectively. The quantity $`b_n`$ is a function of $`n`$, and is evaluated as a root of an equation involving the incomplete gamma function (Wadadekar et al., 1999). It is well approximated by

$$b_n=0.87n0.14.$$

For $`n=4`$, which corresponds to the de Vaucouleurs law, $`b_4=3.33`$. We model the disk intensity by the Freeman (1970) exponential distribution,

$`I_{\mathrm{disk}}(x,y)`$ $`=`$ $`I_d(0)\mathrm{exp}\left(r_{\mathrm{disk}}/r_d\right);`$ (2)
$`r_{\mathrm{disk}}`$ $`=`$ $`\sqrt{x^2+y^2/(1ϵ_d)^2}.`$
The disk is assumed to be intrinsically circular, with the observed ellipticity $`ϵ_d`$ of the disk in the observed image being due to projection effects alone.

The free parameters of the fit are:

1. The central bulge intensity $`I_b(0)`$,
2. The half-light radius of the bulge $`r_e`$,
3. The ellipticity of the bulge $`ϵ_b`$,
4. The bulge shape parameter $`n`$ (Sersic profile),
5. The central intensity of the disk $`I_d(0)`$,
6. The scale length of the disk $`r_d`$, and
7. The ellipticity of the disk $`ϵ_d`$.

Additional parameters such as sky background and position angle are estimated in a preliminary fit and are then held fixed during the final fitting process.

## 3 Modelling the two-dimensional surface brightness of galaxies

We were able to obtain satisfactory fits with a reduced $`\chi _\nu ^{\mathrm{\hspace{0.17em}2}}2`$ for 471 galaxies from our initial sample of 764 galaxies. Any fit with a reduced $`\chi _\nu ^{\mathrm{\hspace{0.17em}2}}>2`$ is considered to be unacceptable and the galaxy excluded from the sample. Such unacceptable fits result mainly from irregularities in the brightness distribution, or from the presence of multiple nuclei or close companions. Such features cannot be represented by the simple and smooth models used in this study.

### 3.1 Morphological classification

The galaxies with satisfactory fits to the surface brightness model are assigned morphological types based on their bulge-to-total luminosity ratio, $`B/T`$. An ideal Elliptical galaxy has $`B/T=1`$ and a late-type disk galaxy has a value of $`B/T`$ close to zero.

We classify a galaxy as early-type when it has $`B/T>0.5`$. This criterion was chosen to be $`B/T>0.4`$ by Balogh et al. (2002), but the conclusions of this study are not sensitive to such a difference in selection. An “Elliptical galaxy” is one with $`B/T0.95`$. A galaxy with $`B/T<0.1`$ is considered as a pure disk galaxy with an exponential luminosity profile, since the bulge in this case is too small to be modelled by the Sersic profile.

### 3.2 Galaxies in X-ray bright, dim and undetected groups

The X-ray luminosities of the parent groups in our sample are taken from Osmond & Ponman (2003), who measured them by fitting two-dimensional $`\beta `$-profiles to ROSAT PSPC observations in the 0.5–2 keV range, extrapolated to estimate the bolometric X-ray luminosity. Point sources were removed from the data before luminosity calculations and the values have then been corrected for flux lost in this process using the surface brightness models.

Since the number of galaxies in each group is small, we obtain statistical results by grouping galaxies according to X-ray brightness of their parent groups. We have divided the groups to be X-ray bright or dim according to their bolometric X-ray luminosity being more or less than the median of the sample, $`L_X=10^{41.7}`$ erg/s. This luminosity refers to the X-ray luminosity of the group plus any central galaxy that might exist. In addition, we have a third category of group (called X-ray undetected) where there is no discernible “group emission”, and all the diffuse emission, if any, can be accounted for as due to individual non-central galaxies in the group.

There are 189, 174 and 108 galaxies (out of 316, 298 and 150) with acceptable profile fits in the categories of X-ray bright, dim and undetected groups respectively. Since the sample covers a wide range of X-ray luminosity, here we seek possible differences in the properties of galaxies that depend on their local environment characterised by $`L_X`$.

### 3.3 Variation in brightness profiles with galaxy environment

Figure 1 shows the distribution of the Sersic parameter for early type galaxies ($`B/T>0.5`$) in X-ray bright, dim and undetected groups. There is no significant difference in the distribution of values of the Sersic parameter $`n`$ among X-ray bright, dim and undetected groups, corroborated by the K-S statistics.

The difference, in the distribution of the half-light radius $`r_e`$ of bulges of early-type disk galaxies, between the various categories of parent groups, is substantial. Figs. 2 and 3 show the distribution of $`r_e`$ and $`M_R`$ in the three sub-samples in this study, indicating that X-ray bright groups contain a larger fraction of galaxies that are bigger and more luminous than their counterparts in X-ray dim groups. This is confirmed by the K-S test performed on the distributions of $`r_e`$ and $`M_R`$. The underdensity of early-types implied by the gap in the histograms for X-ray dim and undetected groups in Fig. 3 around $`M_R=18`$ is due to the dip in the LF of these groups, discussed in further detail by Miles et al. (2003).

The Sersic parameter $`n`$ is often found to be correlated with the absolute magnitude of the Elliptical galaxies (Young & Currie, 1994; Binggeli & Jerjen, 1998), which would indicate that the Sersic parameter is an intrinsic property of a galaxy. Here, we find a highly scattered correlation between $`n`$ and $`M_R`$ for Elliptical galaxies (both E and dE) in groups, with a scatter of about 1.5 mag in $`M_R`$, which is much larger than the scatter (about 0.45 mag in $`M_B`$) found by Young & Currie (1994).

In order to understand the origin of this scatter in the $`nM_R`$ relation, we plot the galaxies as a function of their normalised projected distance from the centre of their parent group (Fig. 4). Since the distribution of the Sersic parameter does not show a significant variation with X-ray properties of the group, we use galaxies belonging to all three sub-samples together to improve statistics. Galaxies in the core of the groups ($`R/R_{\mathrm{vir}}<0.15`$) have higher values for the Sersic parameter (mean value and 95% confidence level of $`1.92_{0.65}^{+0.66}`$), while in the outer regions ($`R/R_{\mathrm{vir}}>0.6`$), galaxies have lower values of $`n`$ and hence are less cuspy (mean value and 95% confidence level of $`1.07_{0.21}^{+0.20}`$). The systematic variation in values of $`n`$ with the distance from the centre of the parent group appears to be independent of how bright the galaxies are.

Trujillo et al. (2002) recently found a triangular distribution between the concentration index of galaxies (related to $`n`$) and the projected two-dimensional local density of galaxies in clusters, which is expected even if $`n`$ (and concentration parameter) and $`M_R`$ were tightly correlated. At higher projected densities, galaxies with both low and high values of concentration index will be binned together due to projection. On the other hand, since galaxies in low-density environments have low $`n`$, at lower projected densities only galaxies with lower values of concentration index will be seen.

Fig. 5 shows a similar triangular distribution between the Sersic parameter $`n`$ and normalised projected distance for the galaxies in our sample. If $`n`$ and $`M_R`$ of Elliptical galaxies were correlated for normal and dwarf Ellipticals, Fig. 5 would seem to indicate that the dwarf Ellipticals are preferentially removed from the core of each group, but enhanced in the outer regions. However, from Fig 4 it appears that the situation is more complicated. We find faint dEs ($`M_R>15.0`$), with relatively higher values of $`n`$, at smaller normalised projected distances from the centre of the groups. This might suggest that these high-$`n`$ dEs are at the core of their groups. This strongly suggests that environmental effects can play a key role in the shape of the light profile (characterised by values of $`n`$), in galaxies and that this parameter is not an intrinsic property of a galaxy.

This also shows that dEs, traditionally known to have near-exponential brightness profiles, can have cuspy cores, and can exist in the central regions of groups. These “nucleated” dEs have been found to be more centrally concentrated in their distribution in detailed studies of nearby clusters, e.g. Virgo and Fornax (Ferguson & Binggeli, 1994). It is interesting to note that such cusps are also known to exist in starburst dEs (e.g. Summers, Stevens, & Strickland, 2001).

## 4 Photometric scaling relations

One of the essential observational ingredients of the study of galaxy evolution is the family of scaling relations such as the Fundamental Plane, characteristic of dynamically relaxed Elliptical galaxies (Djorgovski & Davis, 1987; Dressler et al., 1987). Although the Fundamental Plane for bright Elliptical galaxies is found to be scarcely affected by environmental factors, it may not be so for low mass Ellipticals. While Nieto et al. (1990) showed that dwarf Ellipticals 111Unless otherwise stated, in this paper we will refer to Ellipticals fainter than $`M_R=18`$ as dwarfs (dEs) and those brighter as normal Ellipticals (Es). follow the same trend defined by the Fundamental Plane of bright Elliptical galaxies, albeit with a larger scatter due to structural variation, the dwarf Ellipticals in Virgo cluster studied by Peterson & Caldwell (1993) belong to a different Fundamental Plane.

In the absence of spectroscopic measures of velocity dispersion, which is a key parameter in constructing the Fundamental Plane, one can study those projections of the Fundamental Plane which are purely photometric. In this section, we will study the photometric characteristics of the Elliptical galaxies in our sample (those with $`B/T0.95`$, as defined in §3.1), a subset of our sample of “early-type galaxies” ($`B/T0.5`$)

### 4.1 The Kormendy relation

Elliptical galaxies are known to show a strong correlation between the mean surface brightness $`\mu _e`$, and the half-light radius $`r_e`$, known as the Kormendy relation (KR) $`\mu _e=A_K\mathrm{log}r_e+B_K`$, with the slope $`A_K`$ lying in the range 2.0–3.0 (Kormendy, 1977). The Kormendy relation for cluster Ellipticals shows a r.m.s. scatter of about 0.1 dex, which would result in a 25% error in $`r_e`$ and hence in measuring distance to the galaxies.

Fig. 6 shows the Kormendy relation for Elliptical galaxies in our sample, where galaxies belonging to parent groups which are X-ray bright ($`L_X>10^{41.7}`$ erg/s), X-ray dim and X-ray undetected are plotted with different symbols. The plot shows two distinct populations: the “normal” Ellipticals ($`M_R<18`$, mean KR slope $`2.3\pm 0.4`$) and the dwarf Ellipticals ($`18<M_R<13`$, mean KR slope $`4.1\pm 0.2`$). However, from Fig. 6, it is apparent that instead of a linear KR, we may interpret the distribution of dEs as one of high scatter, truncated by our magnitude limit ($`M_R13`$).

We also note that there are very few galaxies in the zone between the two populations, which results from the significant dip in the luminosity function of groups in the range $`16>M_R>18`$, as seen in Miles et al. (2003).

Comparing with the KR for galaxies in rich clusters, we find that the slope of the KR for our brightest galaxies ($`M_R<19.5`$) 2.75$`\pm `$0.16, compared to 2.43$`\pm `$0.15 for a similar sample in the Coma cluster (Ziegler et al. 1999, plotted as the solid line in Fig. 6).

Extending the KR plot of the Coma cluster to fainter galaxies (e.g. Graham & Guzmán, 2003), Ellipticals in the magnitude range $`15>M_B>18`$ are found to follow a trend, at almost constant $`r_e`$, quite distinct from that of the brighter galaxies, the latter having a slope similar to that found in our sample. Our sample of dwarfs extends to fainter $`M_R13`$ smaller ($`r_e100`$ pc) galaxies, and shows the presence of substantial scatter in both axes. Some contribution to this large scatter may also arise from dependencies on local environment, since we are merging galaxies from 16 groups in a single plot. Our conclusion (in agreement with Graham & Guzmán (2003)) is that the Kormendy relation in itself does not provide strong evidence for any fundamental difference in structure between normal and dwarf Ellipticals.

### 4.2 The Photometric Plane

Recent developments in the study of the structure of galaxies and the modelling of brightness profiles, in particular with the use of the more general Sersic model, provides a unified description of the photometric properties of normal (E) and dwarf Ellipticals (dE), traditionally seen as two distinct families (e.g. Kormendy, 1985). The discontinuity between these two extreme exponential ($`n1`$, dEs) and de Vaucouleurs profiles ($`n4`$, Es) is now filled by both populations possessing a large range of values $`1<n<6`$. The bulges of disk galaxies share the same Photometric Plane with Ellipticals.

The Photometric Plane (PP, Khosroshahi et al., 2000a) is a bivariate relation that links only photometric parameters, obtained by fitting a Sersic model to a galaxy image, $`\mathrm{log}n=A\mathrm{log}r_e+B\mu _0`$. Our simulations show that fitting galaxy brightness distributions with unsuitable models (e.g. the de Vaucouleurs model to a galaxy where $`n4`$), can result in a significant variation in the values of structural parameters. For example, assuming $`n=4`$ in the case of a $`n=2`$ galaxy can result in a change in the value of half-light radius by a factor of three. Thus it is important to allow for a free parameter $`n`$ to account for the cusp and core in Ellipticals or bulges of disk galaxies.

Most scaling relations, e.g. the Kormendy relation or the Fundamental Plane, are insensitive to these variations, since the three parameters of the Sersic model ($`\mu _0`$, $`n`$, $`r_e`$) are correlated and their combinations used in these scaling relations are robust (Bertin et al., 2002). The Photometric Plane, on the other hand, allows us to explore the entire range of these parameters for Elliptical galaxies, and thus can be used more effectively as a tool to examine the connection between structural parameters and the evolutionary history of galaxies.

The Photometric Plane for all Ellipticals in our sample is given by

$`\mathrm{log}n=`$ $`(0.21\pm 0.09)\mathrm{log}r_e`$ (3)
$`(0.074\pm 0.013)\mu _0`$ $`+(1.7\pm 0.3).`$
This equation is similar to that found for the bright Elliptical galaxies in the $`K`$-band study of 42 Coma galaxies (Khosroshahi et al., 2000a). However, the scatter in $`n`$ here is 0.13 $`dex`$, much larger than the scatter found for the Coma cluster sample, which is 0.05 $`dex`$, which can be attributed to the fact that our sample is dominated by dEs. To verify this, we examine the Ellipticals with $`M_R<18`$ in X-ray bright groups (comparable to the Coma sample), to find that the scatter about the Photometric Plane for Es is only 0.06 $`dex`$ in $`\mathrm{log}n`$, compared to 0.10 $`dex`$ for dEs. The fact that the scatter around the the Photometric Plane of normal Es is lower than that of dEs indicates that Es and dEs are structurally different systems even though their $`n`$ values might be similar.

We plot the principal axes ($`K_1`$, $`K_2`$, $`K_3`$) of the Photometric Plane against each other in Fig. 7, where $`K_1=0.99\mathrm{log}n+0.17\mathrm{log}r_e`$, $`K_2=0.16\mathrm{log}n0.92\mathrm{log}r_e0.38\mu _0`$ and $`K_3=0.06\mathrm{log}n+0.37\mathrm{log}r_e0.93\mu _0`$ (directly comparable to the similar plot for Coma Ellipticals in Khosroshahi et al. 2000a). Since this sample has a fainter magnitude limit, it probes the PP for dwarfs more extensively, and demonstrates that the PP for the dEs has a is strikingly low scatter (in contrast, for example, to their scatter in the Kormendy plot). This strongly suggests that dEs might be a more homogeneous class than is normally thought, and follow well-defined structural trends.

To explore a more direct connection with the underlying physical changes in galaxy structure between Es and dEs, we turn to a slightly different version of these correlations in the next section.

### 4.3 The Specific Entropy Plane

Márquez et al. (2001), in their study of cluster Ellipticals, showed that the Photometric Plane is close to a plane of constant specific entropy of galaxies, originally introduced by Lima Neto et al. (1999). Although self-gravitating systems like galaxies do not posess stable states of absolute maximum entropy, the evolution of entropy via two-body relaxation will be a slow process, and Márquez et al. (2000) argue that the violent relaxation process which occurs during galaxy formation may establish a quasi-equilibrium state of constant specific entropy, and they find evidence that both simulated and observed galaxies do indeed show approximately constant entropy per unit mass. In the case of observed galaxies, the properties of the Photometric Plane result from this approximate constancy of specific entropy.

On examination of the surfaces of constant entropy in detail, it is found that they actually exhibit some curvature when represented in the coordinate system of the PP. It can be seen, from the two edge-on views of our plane in Fig. 7, that there are clear indications of curvature in the distribution of the dwarfs. This curvature, rather than stochastic scatter, accounts for much of the larger variance about the plane for dEs (compared to nornal Es) discussed in section 4.2. Could this curvature be indicating that dwarfs really follow a surface of constant specific entropy? In Fig. 8 we plot a representation of the Photometric Plane in terms of $`1/n`$, for our sample, using Eq. (3) above, and overlay a line of constant specific entropy. The latter has been calculated using Eqs. (23) & (24) of Márquez et al. (2001), in conjunction with the relevant relations from their Appendix A. Here, we have only adjusted the value of the specific entropy parameter ($`s_0`$) to find the curve which best matches our data.

The curvature seen in the distribution of dwarfs in Fig. 8 is similar to that the constant $`s_0`$ surface, although the shape does not match in detail. It is also apparent that the normal Es fall above the dwarfs in the Figure, which corresponds to the direction of increasing $`s_0`$. The dashed line corresponds to a specific entropy Plane with $`s_0=25`$ for our sample of Elliptical galaxies dominated by dEs.

## 5 Discussion and Conclusions

This study confirms earlier results that massive Ellipticals show higher values of the Sersic parameter $`n`$, i.e. a cuspier surface brightness profile at their centre (Young & Currie, 1994). From the large scatter in the plot of absolute R magnitude and $`n`$ (Fig. 4), we conclude that the shape of the radial profile, represented by the value of $`n`$, is not uniquely related to the intrinsic luminosity of the galaxy. Furthermore, the dependence of $`n`$ on the location of galaxies within their groups (Fig. 4) shows that environmental dependence is responsible for part of this scatter. Although we find no significant difference between the values of $`n`$ in groups of different X-ray luminosity, we do find that for both dEs and Es (though not the most luminous Es), surface brightness profiles tend to be cuspier (higher $`n`$) for galaxies located in the cores of groups. One explanation for this correlation could be the effects of interaction-induced star formation in triggering nuclear starbursts within more centrally located galaxies.

Normal Ellipticals in groups are found to yield similar Photometric Plane and Kormendy relation parameters to those in rich clusters. There is no difference in the Kormendy relation or Photometric Plane between Ellipticals belonging to the the different classes of X-ray bright, faint or undetected groups. Furthermore, the Kormendy relation is found to be of little value in distinguishing structural differences between dEs and Es, since the KR slope found for dwarfs is dominated by selection effects, most notably the magnitude limit of the sample. The sparse region in the KR plot, separating dEs and Es, is easily interpereted as resulting from the dip seen in the galaxy luminosity function at $`M_R18`$, and discussed in more detail in Miles et al. (2003).

Dwarf and normal Ellipticals are distinguished much more effectively in terms of their location in the Photometric Plane in the space of $`n`$, $`r_e`$ and $`\mu _0`$. The dwarfs inhabit a non-flat surface in this space, with a curvature somewhat akin to a surface of constant specific entropy, as defined by the equations of Márquez et al. (2001), though the match is only approximate. Even if galaxies do have constant specific entropy, a degree of mismatch between observation and the analytical expressions of Márquez et al. (2001) is not surprising, given the numerous assumptions (e.g. mass-traces-light, and isotropic velocities) which have been made in deriving the equations of the isentropic surfaces. The well-defined nature of the locus followed by these dwarfs is a strong indicator of their structural coherence, and the larger scatter that they exhibit about their Photometric Plane is largely due to the curvature of the surface they inhabit, rather than to a large random scatter in their properties. Normal galaxies, on the other hand, are offset relative to the dwarf locus, in a direction corresponding to higher specific entropy. In this context, it is interesting that Márquez et al. (2000) find, from their simulations, that merging tends to increase specific entropy. Hence the offset may indicate that normal Ellipticals, which fall above the dip seen in the galaxy luminosity function, have undergone mergers, whilst dwarfs have not. This supports the idea (Miles et al., 2003) that merging may have played a role in generating the dip in the luminosity function, by converting middle-ranking galaxies into larger ones, whilst dwarfs are largely unaffected due to their smaller merger cross-sections.

## Acknowledgments

We would like to thank Ale Terlevich for help in observation and basic data reduction, John Osmond for making the results of his analysis available for us and Yogesh Wadadekar for the use of his 2D decomposition code.


format: text

text:
# 1 Introduction

## 1 Introduction

The purpose of the present study is twofold. First, we test a novel time-domain algorithm for the Maxwell’s theory of linear, passive (dispersive and absorbing) media. The algorithm is based on (i) the Hamiltonian formalism for evolution differential equations , on (ii) the time leapfrog scheme , and on (iii) the Fourier pseudospectral method  in combination with a change of variables that enhances the spatial grid resolution in designated domains (in the vicinity of medium interfaces) and, thereby, prevents the loss of accuracy due to the aliasing problem of the Fourier transform, while keeping the total spatial grid size fixed , . Boundary conditions at medium interfaces are not fixed in the algorithm, but rather medium parameters are allowed to have spatial discontinuities so that the correct boundary conditions are enforced dynamically , similarly to the wave packet method for quantum mechanical systems with discontinuous potentials. Although, each of these three main ingredients of our algorithm have been used individually in various computational problems in electromagnetism and quantum mechanics, to our knowledge they have never been put together in applications to numerical simulations of propagation of a wideband electromagnetic pulse in passive media and its scattering on targets made of dispersive and absorbing materials. By combining these methods, we have obtained an efficient, true time-domain algorithm that is highly accurate, which is a known virtue of pseudospectral methods of solving partial differential equations .

An essential advantage of time-domain numerical methods is that one can see all the immediate effects the medium and targets have on the propagating wideband wave packet. Yet, a single simulation of scattering of a wideband wave packet is sufficient to determine some basic electromagnetic properties of a target, e.g., transmission and reflection coefficients, in the entire frequency range of the initial wave packet.

It has been reported  that a periodic thin-film metallic grating (either with holes or one-dimensional slits) can transmit more light at certain wavelengths than the projected area of the holes (or slits) in the grating would suggest, while at other wavelengths transmission is almost fully blocked. There is an ongoing discussion about the mechanism of such anomalous transmission. Among suggested mechanisms are formation of dynamical diffraction resonances in periodic metallic structures , surface plasmons whose resonances are enhanced by the array of holes in the grating , and, specific to the slit gratings, open Fabry-Pérot resonant cavities . Since all the approaches produce essentially identical predictions for the transmission and reflection coefficients (the far-field) of the slit grating, it is, perhaps, necessary to have a closer look at the details of the electromagnetic field dynamics in the vicinity of the grating where deviations of theoretical predictions of a particular mechanism from the actual exact solution of Maxwell’s equations might occur (see, e.g., a discussion in a recent work ). This is our second motivation of the present study. We apply our algorithm to scattering of a wideband electromagnetic pulse on various transmission and reflection slit gratings (dispersive metallic and purely dielectric ones). Our time-domain algorithm allows one to observe in details (in real time) formation and decay of long-living resonant excitations of electromagnetic fields as well as formation of resonant transmission and reflection wave fronts. Thanks to the use of the Fourier pseudospectral method in combination with nonuniform grids, an extremely high accuracy of simulations can be achieved in the entire simulation volume and time span at relatively low computational costs. We think that these virtues of our time-domain algorithm would be useful for numerical studies of electromagnetic properties of other nanostructured materials .

The paper is organized as follows. Section 2 is devoted to the Hamiltonian formalism applied to the Maxwell’s theory for general passive linear media. Maxwell’s equations for electromagnetic fields and medium responses are shown to be equivalent to a first order evolution differential equation similar to the Schrödinger equation in quantum mechanics. The wave function is a multidimensional column whose components are electromagnetic fields and auxiliary fields describing the medium response. The Hamiltonian operator is a linear differential operator acting on a Hilbert space spanned by square integrable wave functions. In Section 3, an example of the multi-resonant Lorentz model, which is widely used to describe passive media, is considered in the framework of the Hamiltonian formalism. In Section 4, we establish a relation between the norm of wave functions and electromagnetic energy. The discussion is limited to the Lorentz model and non-dispersive dielectric media. Section 5 is devoted to a general description of our algorithm. In particular, we show how the action of the Hamiltonian on wave functions is defined on a grid by means of the fast Fourier transform method. We prove that the Gauss law is enforced in the grid representation at no extra computational costs in our algorithm. Then we discuss how the sampling efficiency of the fast Fourier method can be enhanced in designated spatial regions (typically at medium interfaces) by changing variables on the grid. The time evolution is done by a modified leapfrog method applied to the Schrödinger equation. We show that the conventional leapfrog method leads to an unstable algorithm for media with absorption and propose a general method to modify the leapfrog scheme to obtain a conditionally stable algorithm. Finally, we give an explicit realization of our algorithm in the case of the multi-resonant Lorentz model. Section 6 contains a detailed description of the actual computational scheme used in our simulations of extraordinary transmission and reflection grating. We analyze the stability of the scheme using our general approach developed in Section 5. Section 7 is devoted to our numerical results. We also compare them with previous theoretical and experimental studies. Section 8 contains a brief conclusion.

## 2 Maxwell’s theory for passive media in the Hamiltonian formalism

Let $`𝐄`$ and $`𝐇`$ be electric and magnetic fields, respectively, $`𝐃`$ and $`𝐁`$ the corresponding inductions, and $`𝐏`$ and $`𝐌`$ the medium polarization and magnetization vectors. Boldface letters denote three-vector fields in $`^3`$. Propagation of an electromagnetic wave packet in passive linear media in the absence of external radiating sources is described by the following set of equations

$`\dot{\psi }^I(t)`$ $`=`$ $`_0\psi ^F(t),`$ (2.1)
$`\psi ^I(t)`$ $`=`$ $`\psi ^F(t)+\psi ^R(t)=\psi ^F(t)+{\displaystyle _0^t}𝑑\tau \chi (t\tau )\psi ^F(\tau ),`$ (2.2)
$`\psi ^I`$ $`=`$ $`\left(\begin{array}{c}𝐃\\ 𝐁\end{array}\right),\psi ^F=\left(\begin{array}{c}𝐄\\ 𝐇\end{array}\right),\psi ^R=\left(\begin{array}{c}𝐏\\ 𝐌\end{array}\right),_0=\left(\begin{array}{cc}0& c\mathbf{}\times \\ c\mathbf{}\times & 0\end{array}\right),`$ (2.3)
where $`c`$ is the speed of light in vacuum, the overdot denotes the partial derivative $`/t`$ with respect to time $`t`$, the spatial argument of the fields, denoted below by $`𝐫`$, is suppressed. For generic anisotropic media, the medium response function $`\chi (t)`$ is regarded as a linear operator (matrix) acting on the components of $`\psi ^F`$. For isotropic media, it is a scalar. The response function is also position dependent for non-homogeneous media. Let $`\psi ^I(0)=\psi ^F(0)`$ be an initial wave packet with finite energy (finite $`𝕃_2(^3)`$ norm). We are interested in a finite norm solution of the initial value problem for Maxwell equations (2.1) subject to the constraint (the Gauss law)

$$\mathbf{}𝐁(t)=\mathbf{}𝐃(t)=0.$$
(2.4)
The response function $`\chi (t)`$ is usually deduced from a microscopic model of the medium in question . Therefore it is natural to assume that $`\chi (t)`$ is a fundamental solution of some linear evolution differential equation so that

$$_t\psi ^R(t)=\omega _p\psi ^F(t),$$
(2.5)
where $`_t`$ is a linear differential operator (a polynomial in $`/t`$) and $`\omega _p`$ describes a coupling between applied electromagnetic fields and matter. In general, $`\omega _g`$ is a position dependent matrix. Causality of the medium response requires that the Fourier transform $`\stackrel{~}{\chi }(\omega )`$ of the response function should have poles only in the lower part of the complex plane of $`\omega `$  so that the Fourier transform of (2.2) reads $`\stackrel{~}{\psi }^R(\omega )=\stackrel{~}{\chi }(\omega )\stackrel{~}{\psi }^F(\omega )`$. By taking the Fourier transform of Eq. (2.5), one finds that $`\stackrel{~}{}(\omega )\stackrel{~}{\psi }^R(\omega )=\omega _p\stackrel{~}{\psi }^F(\omega )`$ and, hence, $`_t=\omega _p^1[\stackrel{~}{\chi }(i/t)]^1`$. If the response function is known from measurements in the frequency domain, components of $`[\stackrel{~}{\chi }(\omega )]^1`$ can always be approximated by a polynomial with sufficient accuracy in a frequency range of interest.

Now the Hamiltonian formalism  can be applied to (2.5) to transform it to an equivalent system of first-order differential equations

$$\dot{\xi }(t)=_M^F\xi (t)+𝒱_{MF}\psi ^F(t),$$
(2.6)
where $`\xi `$ is a column of auxiliary fields which are linear combinations of the response field and its time derivatives, save for the one of the highest order. There exists a linear operator $``$ such that $`\psi ^R(t)=\xi (t)`$. Its explicit form depends on the details of going over to the Hamiltonian formalism. One can, for instance, identify the first component of $`\xi `$ with $`\psi ^R`$. In this case, $``$ projects the column $`\xi `$ onto its first component. The response function can be expressed through the operators $`_M^F`$ and $`𝒱_{MF}`$ and $``$ by making use of the fundamental solution of Eq. (2.6)

$$\chi (t)=\theta (t)e^{_M^Ft}𝒱_{MF},$$
(2.7)
where $`\theta (t)`$ is the Heaviside function. Equation (2.7) can be regarded as a condition on possible choices of the operators $`_M^F`$ and $`𝒱_{MF}`$ and $``$.

Applying $``$ to (2.6) we find

$$\dot{\psi }^R(t)=_M^F\xi (t),$$
(2.8)
where the relation $`𝒱_{MF}=0`$ has been used. It is not hard to be convinced that the latter relation holds when the first component of $`\xi `$ coincides with $`\psi ^R`$. Any other choice of $`\xi `$ can be obtained by a canonical transformation  which is a linear nonsingular transformation of auxiliary fields,

$$\xi 𝒮_M\xi ,det𝒮_M0.$$
(2.9)
According to (2.6) and (2.7), $`𝒮_M^1`$, $`_M^F𝒮_M_M^F𝒮_M^1`$, and $`𝒱_{MF}𝒮_M𝒱_{MF}`$ and, therefore, the condition $`𝒱_{MF}=0`$ remains true in the new basis of auxiliary fields. Denoting $`𝒱_{FM}=_M^F`$ and substituting (2.2) and (2.8) into (2.1) Maxwell’s evolution equations (2.1) can be written in the Hamiltonian form

$$\dot{\psi }^F(t)=_0\psi ^F(t)+𝒱_{FM}\xi (t).$$
(2.10)
Finally, the electromagnetic and auxiliary fields are unified into one column (wave function) so that Eqs. (2.10) and (2.6) can be represented as a single first-order evolution equation

$`\dot{\mathrm{\Psi }}^F(t)`$ $`=`$ $`^F\mathrm{\Psi }^F(t),`$ (2.11)
$`\mathrm{\Psi }^F`$ $`=`$ $`\left(\begin{array}{c}\psi ^F\\ \xi \end{array}\right),^F=\left(\begin{array}{cc}_0& 𝒱_{FM}\\ 𝒱_{MF}& _M^F\end{array}\right).`$ (2.12)
The index $`F`$ indicates that electromagnetic fields are used as independent electromagnetic degrees of freedom. We shall refer to (2.12) as to a field representation. Accordingly, an induction representation is obtained by the similarity transformation

$$\mathrm{\Psi }^I=\left(\begin{array}{c}\psi ^I\\ \xi \end{array}\right)=𝒮\mathrm{\Psi }^F,𝒮=\left(\begin{array}{cc}1& \\ 0& 1\end{array}\right),^I=𝒮^F𝒮^1.$$
(2.13)
The corresponding blocks of $`^I`$ have the form

$`_I`$ $`=`$ $`_0,𝒱_{MI}=𝒱_{MF},`$ (2.14)
$`𝒱_{IM}`$ $`=`$ $`𝒱_{FM}+_M^F_0=_0,`$ (2.15)
$`_M^I`$ $`=`$ $`_M^F𝒱_{MF}.`$ (2.16)
In what follows we denote wave functions by $`\mathrm{\Psi }^Q(t)`$ with $`Q`$ being the representation index, $`F`$ or $`I`$. This completes construction of the Hamiltonian representation of Maxwell’s theory for linear passive media.

Boundary conditions at medium interfaces and possible scatterers (targets) are not imposed on electromagnetic fields, but rather they are enforced dynamically by allowing medium parameters to be discontinuous functions. The fundamental solution of (2.11)

$$\mathrm{\Psi }^F(t)=e^{t^F}\mathrm{\Psi }^F(0),t0,$$
(2.17)
is well defined for discontinuous “potentials” $`𝒱_{MF}`$ and $`𝒱_{FM}`$, for example, by means of the Kato-Trotter product formula used in the path integral representation of (2.17) as shown in .

## 3 The Lorentz model

The Hamiltonian formalism for the Lorentz model has been used in  to develop a finite differencing algorithm to study an electromagnetic pulse propagation in Lorentz media. Here we derive an explicit form of the Hamiltonian for the Lorentz model which is used in Section 5 for the stability analysis of our algorithm. The Lorentz model of a passive medium is based on the assumption that the medium magnetization is zero, $`𝐌=0`$, while the medium polarization is described by a set of decoupled second-order differential equations

$$\ddot{𝐏}_a+2\gamma _a\dot{𝐏}_a+\omega _a^2𝐏_a=\omega _{pa}^2𝐄,𝐏=\underset{a=1}{\overset{N}{}}𝐏_a.$$
(3.1)
where $`\omega _a`$ are resonant frequencies, $`\gamma _a`$ are damping coefficients, and $`\omega _{pa}`$ are plasma frequencies. As has been pointed out, no boundary conditions are imposed on electromagnetic fields at medium and/or target interfaces. Instead, the coupling constants $`\omega _{pa}=\omega _{pa}(𝐫)`$ are allowed to have discontinuities at medium interfaces, or, from the physical point of view, they remain smooth but change rapidly, $`\lambda _w|\mathbf{}\omega _p|/\omega _p>>1`$, at the interface, where $`\lambda _w`$ is a typical wave length of the incident wave packet. The initial value problem is solved in the space of square integrable wave functions. Initial conditions for the response field are $`𝐏(t=0)=\dot{𝐏}(t=0)=0`$.

Consider $`2N`$ real vector fields, $`𝝃_j`$, $`j=1,2,\mathrm{},2N`$, such that

$`𝐏_a`$ $`=`$ $`(\omega _{pa}/\omega _a)𝝃_{2a1},`$ (3.2)
$`\dot{𝝃}_{2a1}`$ $`=`$ $`\omega _a𝝃_{2a},\dot{𝝃}_{2a}=2\gamma _a𝝃_{2a}\omega _a𝝃_{2a1}+\omega _{pa}𝐄.`$ (3.3)
Thus, the original system of second order equations has been converted into the first order system. The operator $``$ is defined by (3.2). After simple algebraic transformations, we infer

$`𝒱_{FM}`$ $`=`$ $`(𝒱_{FM1},𝒱_{FM2},\mathrm{},𝒱_{FMN}),𝒱_{FMa}=\left(\begin{array}{cc}0& \omega _{pa}\\ 0& 0\end{array}\right),`$ (3.4)
$`𝒱_{MF}`$ $`=`$ $`𝒱_{FM}^{},`$ (3.5)
$`_M^F`$ $`=`$ $`\mathrm{𝖽𝗂𝖺𝗀}(_{M1}^F,_{M2}^F,\mathrm{},_{MN}^F),_{Ma}^F=\left(\begin{array}{cc}0& \omega _a\\ \omega _a& 2\gamma _a\end{array}\right),`$ (3.6)
where $`\mathrm{𝖽𝗂𝖺𝗀}`$ indicates that the corresponding matrix is block-diagonal with blocks listed in the order from the upper left to lower right corners. Note that the matrices $`𝒱_{FMa}`$ and $`_{FMa}`$ act on a six-dimensional column $`\xi _a`$ composed of two vectors $`𝝃_{2a1}`$ and $`𝝃_{2a}`$. Therefore they should be understood as composed of $`3\times 3`$ blocks. Each block is obtained by multiplying the unit matrix by the number indicated in place of the block in (3.4) and (3.6).

Another convenient way to introduce the Hamiltonian formalism is to use $`N`$ complex vector fields $`𝜻_a`$ which satisfy the first order differential equation

$$\dot{𝜻}_a=\lambda _a𝜻_ai\omega _{pa}𝐄,𝐏_a=\frac{\omega _{pa}}{2\nu _a}\left(𝜻_a+\overline{𝜻}_a\right),$$
(3.7)
where $`\lambda _a=\gamma _a+i\nu _a`$ and $`\nu _a=\sqrt{\omega _a^2\gamma _a^2}`$. This representation is defined only if $`\gamma _a<\omega _a`$ (i.e., the attenuation is not high). From the numerical point of view, solving a decoupled system of $`N`$ first order differential equation and taking complex conjugation (denoted here by an over bar) is less expensive than solving the original system of second-order differential equations for the medium polarization. In general, there is always a freedom of choosing a new basis for the auxiliary field space (2.9). If the evolution operator $`\mathrm{exp}(t^Q)`$ is computed in one basis, it can be computed in another basis by a suitable similarity transformation. This is an important observation because the auxiliary field basis can be chosen in a way that facilitates computation of the evolution operator, e.g., to speed up simulations. For instance, in the complex representation (3.7), the matter Hamiltonian $`_M^F`$ is diagonal. The corresponding transformation of auxiliary fields is given by

$$\left(\begin{array}{c}𝝃_{2a1}\\ 𝝃_{2a}\end{array}\right)=\frac{1}{2\nu _a}\left(\begin{array}{cc}\omega _a& \omega _a\\ \lambda _a& \overline{\lambda }_a\end{array}\right)\left(\begin{array}{c}𝜻_a\\ \overline{𝜻}_a\end{array}\right)𝒮_M\left(\begin{array}{c}𝜻_a\\ \overline{𝜻}_a\end{array}\right).$$
(3.8)
To transform the whole system into this representation, the Hamiltonian $`^F`$ is replaced by $`𝒮^1^F𝒮`$ and the wave function $`\mathrm{\Psi }^F`$ by $`𝒮\mathrm{\Psi }^F`$ where $`𝒮`$ is block-diagonal with the unit matrix in the upper left (field) corner and with $`𝒮_M`$ in the lower right (matter) corner.

## 4 Energy and the norm of state vectors

Accuracy and convergence of a numerical algorithm is defined relative to some norm. Let us discuss the choice of a norm in the space spanned by wave functions $`\mathrm{\Psi }^Q(t)`$. The discussion is limited to the Lorentz model and the case of a nonhomogeneous, nondispersive medium (dielectric) which are used in our numerical simulations.

### I The Lorentz model

Consider a multi-resonant Lorentz model with no attenuation $`\gamma _a=0`$. The field and matter evolution equations can be obtained from the variational principle for the action

$$S=𝑑tL=𝑑t𝑑𝐫\left[\frac{1}{2}\left(𝐄^2𝐁^2\right)+\frac{1}{2}\underset{a}{}\left(\dot{\mathit{\vartheta }}_a^2\omega _a^2\mathit{\vartheta }_a^2\right)+𝐏𝐄\right],$$
(4.1)
where the polarization of the medium is expressed via matter fields as $`𝐏=_a\omega _{pa}\mathit{\vartheta }_a`$. The electromagnetic degrees of freedom are described by the vector and scalar potentials, respectively, $`𝐀`$ and $`\phi `$. The fields are defined by $`𝐄=\mathbf{}\phi \dot{𝐀}`$ and $`𝐁=\mathbf{}\times 𝐀`$. The units are chosen in this Section so that $`c=1`$. The least action principle for the scalar potential $`\phi `$ leads to the Gauss law, $`\mathbf{}𝐃=0`$, for the vector potential $`𝐀`$ to the Maxwell’s equation, $`\dot{𝐃}=\mathbf{}\times 𝐁`$, and for the matter fields $`\mathit{\vartheta }_a`$ to the medium polarization evolution equation (3.1) of the Lorentz model with no attenuation, $`\gamma _a=0`$. The second Maxwell’s equation and the Gauss law for the magnetic field follows from the relation $`𝐁=\mathbf{}\times 𝐀`$ by taking its time derivative and divergence, respectively.

The energy of the system coincides with the canonical Hamiltonian which is obtained by the Legendre transformation  of the Lagrangian $`L`$ for the velocities $`\dot{𝐀}`$ and $`\dot{\mathit{\vartheta }}_a`$. The canonical momenta are $`𝝅_a=\delta L/\delta \dot{\mathit{\vartheta }}_a=\dot{\mathit{\vartheta }}_a`$ and $`𝚷=\delta L/\delta \dot{𝐀}=𝐄𝐏=𝐃`$. Doing the Legendre transformation, we find the canonical Hamiltonian (energy) of the system

$$E(t)=𝑑𝐫\left(\underset{a}{}𝝅_a\dot{\mathit{\vartheta }}_a+𝚷\dot{𝐀}\right)L=\frac{1}{2}𝑑𝐫\left[𝐄^2+𝐁^2+\underset{a}{}\left(𝝅_a^2+\omega _a^2\mathit{\vartheta }_a^2\right)\right],$$
(4.2)
where the Gauss law $`\mathbf{}𝐃=0`$ has been used. The energy conservation follows directly from the Noether theorem  applied to the time translation symmetry of the action (4.1), $`\dot{E}(t)=0`$. Equation (4.2) becomes the conventional expression for the electromagnetic energy in a passive medium  when $`𝝅_a`$ and $`\mathit{\vartheta }_a`$ are replaced by the corresponding solutions of the equations of motion with initial conditions $`𝝅_a(t=0)=\mathit{\vartheta }_a(t=0)=0`$.

An important observation is that the Noether integral of motion (4.2) coincides with the norm squared of the corresponding state vector

$$E=\frac{1}{2}𝑑𝐫\mathrm{\Psi }^F\mathrm{\Psi }^F(\mathrm{\Psi }^F,\mathrm{\Psi }^F)=(\mathrm{\Psi }^I,\mu \mathrm{\Psi }^I),$$
(4.3)
where $`\mu =𝒮^1𝒮^1`$. This follows from the fact that, if we identify $`𝝃_{2a}=𝝅_a`$ and $`𝝃_{2a1}=\omega _a\mathit{\vartheta }_a`$, the Hamiltonian equations of motion, $`\dot{\mathit{\vartheta }}_a=\delta E/\delta 𝝅_a`$ and $`\dot{𝝅}_a=\delta E/\delta \mathit{\vartheta }_a`$, coincide with (3.3) when $`\gamma _a=0`$. Note that canonically conjugated electromagnetic variables are $`𝐀`$ and $`𝐃`$. Therefore, the coupling between the electromagnetic and matter degrees of freedom in the Hamiltonian equations of motion is generated by the term $`𝐄^2=(𝐃𝐏)^2`$ in (4.2). Thus, in the absence of attenuation, the norm of the state vector is proportional to the wave packet electromagnetic energy and, hence, is conserved.

The norm conservation also follows from anti-Hermiticity of the Hamiltonian $`^F=^F`$ if $`\gamma _a=0`$, while (4.3) establishes a relation between the electromagnetic energy and the norm. In the induction representation, the norm in the measure space, defined by the operator $`\mu `$ in (4.3), is also conserved by construction. Consequently, the Hamiltonian is anti-Hermitian relative to the measure space scalar product, $`^I\mu =\mu ^I`$.

The norm (energy) conservation can be used to control numerical convergence, especially when the aliasing problem in the fast Fourier transform is present, e.g., when parameters of the medium are discontinuous functions in space. In a properly designed algorithm the loss of energy (norm) due to attenuation should be controlled by the symmetric part of the Hamiltonian operator

$$\dot{E}(t)=\underset{a}{}\gamma _a𝑑𝐫𝝃_{2a}^2(t)\frac{1}{2}(\mathrm{\Psi }^F(t),𝒱_\gamma ^F\mathrm{\Psi }^F(t))0,$$
(4.4)
where $`𝒱_\gamma ^F=𝒱_\gamma ^F=(^F+^F)/20`$ (a negative semidefinite operator) which is, in this case, a diagonal matrix with nonpositive elements.

### II Nondispersive media

If the medium in question does not have dispersion and absorption, the formalism is simplified. Let $`\epsilon =\epsilon (𝐱)`$ be the dielectric constant of the medium, $`𝐃=\epsilon 𝐄`$. If the medium is not isotropic, then $`\epsilon `$ is symmetric positive definite $`3\times 3`$ matrix everywhere in space. We rewrite Maxwell’s equations in the form

$$\dot{\psi }^I(t)=_G\psi ^I(t),_G=\left(\begin{array}{cc}0& c\mathbf{}\times \\ c\mathbf{}\times (\epsilon ^1)& 0\end{array}\right),$$
(4.5)
where the parenthesis in $`(\epsilon ^1)`$ mean that the induction is first multiplied by $`\epsilon ^1`$ and then the curl of the resulting vector field is computed. Consider the scalar product

$$(\psi _1^I,\mu _\epsilon \psi _2^I)=𝑑𝐫\psi _1^I\mu _\epsilon \psi _2^I,\mu _\epsilon =\left(\begin{array}{cc}\epsilon ^1& 0\\ 0& 1\end{array}\right).$$
(4.6)
The Hamiltonian is anti-Hermitian with respect to this measure space scalar product, $`_G^{}\mu _\epsilon =\mu _\epsilon _G`$. Therefore the corresponding norm is preserved in the time evolution generated by $`\mathrm{exp}(t_G)`$, that is, $`(\psi ^I(t),\mu _\epsilon \psi ^I(t))=(\psi ^I(0),\mu _\epsilon \psi ^I(0))`$. The electromagnetic energy of the wave packet is conserved because it is proportional to the measure space norm of the initial state vector.

## 5 The algorithm

### I The grid representation

Consider an equidistantly spaced finite grid with periodic boundary conditions. Let $`\mathrm{\Delta }r`$ be the grid step and $`𝐧`$ be a vector with integer valued components. Then the dynamical variables are projected onto the grid by taking their values at grid points $`𝐫=𝐧\mathrm{\Delta }r`$, that is, the wave function is replaced by a vector (column) whose dimension is determined by the grid size, $`\mathrm{\Psi }^Q(𝐫)\mathrm{\Psi }^Q(𝐧\mathrm{\Delta }r)\mathrm{\Psi }_𝐧^Q`$. A cubic grid is assumed. Consider a discrete Fourier transformation associated with the grid, $`\stackrel{~}{\mathrm{\Psi }}^Q(𝐧k_0)=_𝐧^{}_{\mathrm{𝐧𝐧}^{}}\mathrm{\Psi }_𝐧^{}^Q\stackrel{~}{\mathrm{\Psi }}_𝐧^Q`$, where $`^{}=^{}=1`$. The reciprocal lattice step is $`k_0=2\pi /\mathrm{\Delta }r`$. The grid spatial size $`L`$ and step $`\mathrm{\Delta }r`$ must be chosen so that the Fourier transform of the initial wave packet has support within the region $`k[k_{min},k_{max}]`$ where $`k=|𝐤|`$, $`k_{max}=k_0`$ and $`k_{min}=2\pi /L`$. The action is of any position dependent operator $`𝒱^Q=𝒱^Q(𝐫)`$ is defined by

$$𝒱^Q(𝐫)\mathrm{\Psi }_t^Q(𝐫)𝒱^Q(𝐧\mathrm{\Delta }r)\mathrm{\Psi }_t^Q(𝐧\mathrm{\Delta }r)\underset{𝐧^{}}{}𝒱_{\mathrm{𝐧𝐧}^{}}^Q\mathrm{\Psi }_𝐧^{}^Q,$$
(5.1)
where $`𝒱_{\mathrm{𝐧𝐧}^{}}^Q=\delta _{\mathrm{𝐧𝐧}^{}}𝒱^Q(𝐧\mathrm{\Delta }r)`$ is a diagonal matrix. Let $`_0^Q=_0^Q(\mathbf{})`$ depend only on the $`\mathbf{}`$ operator. The projection of its action onto the grid is then defined via the discrete Fourier transform

$$_0^Q(\mathbf{})\mathrm{\Psi }_t^Q(𝐫)|_{𝐫=𝐧\mathrm{\Delta }r}\underset{𝐧^{}}{}_{0\mathrm{𝐧𝐧}^{}}^Q\mathrm{\Psi }_𝐧^{}^Q,_{0\mathrm{𝐧𝐧}^{}}^Q=\underset{𝐧^{\prime \prime }}{}\left(^{}\right)_{\mathrm{𝐧𝐧}^{\prime \prime }}_0^Q(i𝐧^{\prime \prime }k_0)\left(\right)_{𝐧^{\prime \prime }𝐧^{}}.$$
(5.2)
The projection (5.2) is performed by the fast Fourier transform method.

In what follows, the action of a product of $`𝒱^Q`$ and $`_0^Q`$ on any state vector is understood as multiplication of $`\mathrm{\Psi }_𝐧^Q`$ by the corresponding matrices, defined in (5.1) and (5.2), in the order specified in the product. The main advantage of using the Fourier basis is the exponential convergence (versus the polynomial one in finite differencing schemes)  as the grid size increases, which allows one to substantially increase the accuracy of simulations.

### II The Gauss law

Another advantage of the Fourier basis is that the Gauss law is enforced at no extra computational cost. In the grid representation defined above, the Gauss law (2.4) requires that the Fourier transforms of the inductions $`\stackrel{~}{𝐃}(𝐤)`$ and $`\stackrel{~}{𝐁}(𝐤)`$ remains perpendicular to the reciprocal grid vector $`𝐤=𝐧k_0`$ at any moment of time. In our algorithm, as we shall show shortly, the time evolution is generated by applying powers of the Hamiltonian to the wave function. In the induction representation, the action of powers of $`^I`$ always produces the cross product $`𝐤\times \stackrel{~}{𝐂}(𝐤)`$, for some $`\stackrel{~}{𝐂}(𝐤)`$ regular at $`𝐤=0`$, in the entries of $`\mathrm{\Psi }^I`$ that correspond to the electromagnetic inductions. Hence, in the grid representation the wave function $`(^I)^m\mathrm{\Psi }^I`$ satisfies the Gauss law for any power $`m`$ because of the trivial identity $`𝐤(𝐤\times \stackrel{~}{𝐂}(𝐤))=0`$ valid for any vector $`𝐤`$ of the reciprocal grid.

It should be noted that a high accuracy of the Gauss law is essential to achieve a high accuracy of simulated electromagnetic fields near medium interfaces.

### III Improving sampling efficiency by changing variables

As is well known from the Fourier analysis, the convergence rate can be affected for functions which have discontinuities . The latter is, unfortunately, the case in electromagnetic scattering problems . Suppose there is an interface between two media. It can be deduced from the dynamical Maxwell’s equations that the components of electric and magnetic fields, $`𝐄`$ and $`𝐇`$, tangential to the interface must be continuous, provided there is no surface electric current on the interface. From the Gauss law it follows that components of the inductions, $`𝐃`$ and $`𝐁`$, normal to the interface must be continuous, provided there is no surface charge on the interface. In contrast, normal components of the fields and tangential components of the inductions can be discontinuous. Their discontinuities are determined by discontinuities of medium parameters (e.g., discontinuities in plasma frequencies in Lorentz models). Therefore, in either the induction or field representation, there are components which suffer discontinuities at the interface. The only way to cope with the problem, while keeping the use of the Fourier basis, is to make the grid finer . This would lead to a substantial waste of computational resources because the conventional fast Fourier transform is defined on a uniform periodic grid, while the sampling efficiency should only be enhanced in the neighborhood of medium interfaces. The use of wavelet bases might be helpful for such a task  in time domain algorithms. Here we retain the Fourier basis, and increase the sampling efficiency by a change of variables.

The basic idea can be understood with a one-dimensional example. Let $`z`$ be a physical coordinate. Consider a change of variables defined by $`z=f(y)`$ where $`y`$ is an auxiliary coordinate. An equidistant grid $`y_n=n\mathrm{\Delta }y`$, with $`n`$ being integers, of the auxiliary coordinate generates a non-uniform grid of the the physical coordinate, $`z_n=f(n\mathrm{\Delta }y)`$. Assuming $`\mathrm{\Delta }y`$ to be sufficiently small and $`f(y)`$ sufficiently smooth, the physical grid spacing can be approximated as

$$\mathrm{\Delta }z_n=z_{n+1}z_n\mathrm{\Delta }yf^{}(n\mathrm{\Delta }y).$$

So, if $`f^{}(y)=1`$, then $`\mathrm{\Delta }z_n=\mathrm{\Delta }y`$ and the grid is equidistant. By making the derivative $`0<f^{}(y)1`$ in some designated areas, one can achieve a desired local grid density in the physical space, while keeping the total grid size fixed. For example, if it is necessary to increase the sampling efficiency in the vicinity $`z=0`$, one can take $`f^{}(y)=1a_0[1+b_0^2y^2]^1`$, where $`0<a_0<1`$, and , hence,

$$f(y)=y\frac{a}{b}\mathrm{tan}^1(by)yg(y,a,b).$$

By adjusting parameters $`a`$ and $`b`$, the local grid density can be changed as desired. Consequently, if the sampling efficiency is to be enhanced at several points $`y_i`$, a suitable change of variables can be of the form $`f(y)=y_ig(yy_i,a_i,b_i)`$. The fast Fourier transform is applied on a uniform grid of the auxiliary variable.

Upon the change of variables, the integration measure in the scalar product and the derivative transform, respectively, $`dz=dyf^{}(y)`$ and $`_z=[f^{}(y)]^1_y`$. When projected on a uniform grid of the physical coordinate by means of (5.2), the derivative operator $`_z`$ becomes an anti-Hermitian matrix. When the rule (5.2) is applied on a uniform grid of the auxiliary coordinate $`y`$, the derivative operator $`_z`$ is no longer represented by an anti-Hermitian matrix, although it is still an anti-Hermitian linear operator, but in the measure space where the scalar product is defined with the weight $`f^{}(n\mathrm{\Delta }y)`$ at each lattice cite. From the numerical point of view it is convenient to have an explicitly anti-Hermitian matrix representation of the derivative on the grid. Due to round-off errors, the exact anti-Hermiticity of $`_z`$ in the measure space can be violated in the grid representation, which, in turn, may lead to numerical instabilities of simulations. Note that if a computed matrix is known to be anti-Hermitian, then only half of its elements is to be computed, while the other elements are restored by the symmetry. In the measure space, an explicit anti-Hermiticity of a linear operator is much more difficult to maintain in numerical simulations because the symmetry relation between matrix elements depends on the scalar product measure. For this reason, the wave function is rescaled by the square root of the Jacobian , $`\mathrm{\Psi }\sqrt{f^{}}\mathrm{\Psi }`$ so that the integration measure becomes $`dy`$ leading to the conventional Euclidean scalar product in the grid representation (with the uniform unit weight at each grid site). In such a representation, the derivative operator $`_z[f^{}]^{1/2}_y[f^{}]^{1/2}`$ becomes again an explicitly anti-Hermitian matrix in the grid representation defined by the rules (5.1) and (5.2).

### IV A modified temporal leapfrog scheme

A conventional temporal leapfrog method of solving the Schrödinger equation is based on the iteration scheme

$$\mathrm{\Psi }^Q(t+\mathrm{\Delta }t)=\mathrm{\Psi }^Q(t\mathrm{\Delta }t)+2\mathrm{\Delta }t^Q\mathrm{\Psi }^Q(t),$$
(5.3)
so that the wave function in the consecutive time moment is computed from the wave function at two previous moments of time, where $`\mathrm{\Delta }t`$ is the time step. The action of the Hamiltonian $`^Q`$ is computed by pseudospectral methods, in particular, by means of the Fourier basis and the fast Fourier transform in our algorithm. The leapfrog algorithm is conditionally stable for an anti-Hermitian $`^Q`$, which is true for nondispersive media, but not the case for media with absorption.

To investigate the stability, let us introduce the amplification matrix for the leapfrog algorithm, $`\mathrm{\Psi }(t+\mathrm{\Delta }t)=𝒢(\mathrm{\Delta }t)\mathrm{\Psi }(t)`$ (the representation index, $`Q`$, is suppressed for a moment). From (5.3) it follows that $`𝒢(\mathrm{\Delta }t)`$ satisfies the equation $`𝒢^2(\mathrm{\Delta }t)2\mathrm{\Delta }t𝒢(\mathrm{\Delta }t)1=0`$ which has two solutions

$$𝒢^{(\pm )}(\mathrm{\Delta }t)=\mathrm{\Delta }t\pm \sqrt{1+^2\mathrm{\Delta }t^2}.$$
(5.4)
The stability of the algorithm requires that the energy norm of both approximate solutions $`[𝒢^{(\pm )}(\mathrm{\Delta }t)]^n\mathrm{\Psi }(0)`$ must be uniformly bounded for $`n>0`$. The necessary condition (but not sufficient) is the von Neumann condition that the spectral radius of the amplification matrix does not exceed 1. If a complex number $`Re^{i\phi }`$ is an eigenvalue of $`\mathrm{\Delta }t`$, then the von Neumann condition implies that $`\phi =\pm \pi /2`$ and $`R^21`$. In other words, eigenvalues of $`\mathrm{\Delta }t`$ must be imaginary and their magnitude should not exceed 1. If in addition we demand that the Hamiltonian is diagonalizable, then a conditional stability can be achieved for sufficiently small $`\mathrm{\Delta }t`$. Indeed, in this case there exists an non-singular $`𝒮`$ such that

$$=𝒮^1_S𝒮,_S^{}=_S,$$
(5.5)
and $`\mathrm{\Delta }t_S`$ satisfies the von Neumann condition. Since $``$ and $`_S`$ have the same eigenvalues, the amplification matrices for $``$ and $`_S`$ are related by the same similarity transformation (5.5). Hence, the norm of a wave function obtained by the action of powers of (5.4) on an initial wave function is uniformly bounded. Note also that the Hamiltonian (5.5) is anti-Hermitian relative to the measure space scalar product, $`^{}\mu =\mu `$, where $`\mu =𝒮^1𝒮^1`$. The stability can also be proved via the equivalence of the conventional Euclidean norm and the $`\mu `$-norm .

In the case of nondispersive media, the Hamiltonian $`_G`$ is anti-Hermitian in the measure space, and the von Neumann condition is fulfilled if

$$\mathrm{\Delta }tck_{max}^\epsilon 1,$$
(5.6)
where $`k_{max}^\epsilon `$ the maximal norm of all wave vectors in the medium which can be estimated by $`\sqrt{\rho (\epsilon )}k_{max}`$ with $`k_{max}`$ being the maximal wave vector of the initial pulse in vacuum and $`\rho (\epsilon )`$ the maximal spectral radius of the symmetric matrix $`\epsilon (𝐱)`$ over $`𝐱`$ (or simply the maximum of $`\epsilon (𝐱)`$ if the medium is isotropic). This can be understood from the following principle . A finite difference scheme with variable coefficients is stable if all the corresponding schemes with frozen (i.e., fixed to a particular value everywhere in space) coefficients are stable.

The von Neumann condition cannot be met if absorption is present because eigenvalues of the Hamiltonian must have real parts in order to account for exponential attenuation of field amplitudes. To circumvent this difficulty, the leapfrog scheme is modified in the following way . We assume the Hamiltonian to be diagonalizable. The lack of eigenvectors of the Hamiltonian typically leads to solutions whose amplitudes grow polynomially in time . This feature cannot be present in a physically reasonable model of passive media. So our assumption is justified from the physical point of view and, yet, the Lorentz model Hamiltonian is indeed diagonalizable. Let

$$=_0+𝒱,(\mathrm{\Psi },𝒱\mathrm{\Psi })0$$
(5.7)
for any $`\mathrm{\Psi }`$ in the Hilbert space and the von Neumann condition is satisfied for $`_0`$, i.e., $`_0`$ has imaginary eigenvalues. The split (5.7) can be achieved in many ways. For instance, $`_0`$ can be obtained from $``$ by setting all parameters responsible for attenuation to zero. In a physically acceptable model, the energy must be conservative in an absorption free medium and so must be the energy norm of wave vectors, and, therefore, the corresponding Hamiltonian must be anti-Hermitian (relative to the energy induced scalar product). In the Lorentz model this is easily seen in the field representation if we set $`_0=|_{\gamma _a=0}`$ which is explicitly anti-Hermitian (see Section 3). Then $`𝒱`$ is diagonal with matrix elements being zeros and $`2\gamma _a`$. Another possibility is to identify $`_0`$ with the Hamiltonian in the vacuum, then $`𝒱=_0`$ must be negative semidefinite in order to model exponential attenuation in passive media. Finally, one can also split the Hamiltonian into the sum of Hermitian and anti-Hermitian parts.

After choosing a suitable split (5.7) we make a substitution $`\mathrm{\Psi }(t)=\mathrm{exp}(t𝒱)\mathrm{\Phi }(t)`$ in the original evolution equation (2.11). The new wave function $`\mathrm{\Phi }(t)`$ satisfies the equation with a time dependent Hamiltonian

$$\dot{\mathrm{\Phi }}(t)=e^{t𝒱}_0e^{t𝒱}\mathrm{\Phi }(t)(t)\mathrm{\Phi }(t),$$
(5.8)
to be solved with the same initial condition $`\mathrm{\Phi }(0)=\mathrm{\Psi }(0)`$. Applying the leapfrog method to (5.8) we get $`\mathrm{\Phi }(t+\mathrm{\Delta }t)=\mathrm{\Phi }(t\mathrm{\Delta }t)+2\mathrm{\Delta }t(t)\mathrm{\Phi }(t)`$ valid up to $`O(\mathrm{\Delta }t^3)`$. Returning to the initial variables, we arrive at the following recurrence relation

$$\mathrm{\Psi }(t+\mathrm{\Delta }t)=(2\mathrm{\Delta }t)\mathrm{\Psi }(t\mathrm{\Delta }t)+2\mathrm{\Delta }t(\mathrm{\Delta }t)_0\mathrm{\Psi }(t),$$
(5.9)
where $`(\mathrm{\Delta }t)=\mathrm{exp}(\mathrm{\Delta }t𝒱)`$. The amplification matrix, $`\mathrm{\Psi }(t+\mathrm{\Delta }t)=𝒢_{}(\mathrm{\Delta }t)\mathrm{\Psi }(t)`$, for the recurrence (5.9) satisfies the equation

$$𝒢_{}(\mathrm{\Delta }t)=(2\mathrm{\Delta }t)𝒢_{}^1(\mathrm{\Delta }t)+2\mathrm{\Delta }t(\mathrm{\Delta }t)_0.$$
(5.10)
A deviation of the approximate solution $`𝒢_{}^n(\mathrm{\Delta }t)\mathrm{\Psi }(0)`$ from the exact solution relative to the energy norm is of order $`\mathrm{\Delta }t^2`$ for any $`n>0`$. Thus, the scheme is convergent and, hence, a conditional stability exists for a sufficiently small $`\mathrm{\Delta }t>0`$ according to a general theorem of Kantorovich  that establishes a general equivalence between convergence and conditional stability. The conditional stability of (5.9) can be understood from the following observation. Solving (5.10) by perturbation theory in $`\mathrm{\Delta }t`$, it is not hard to find that

$$𝒢_{}(\mathrm{\Delta }t)𝒢_V(\mathrm{\Delta }t)=\mathrm{\Delta }t^3𝒦(\mathrm{\Delta }t),𝒢_V(\mathrm{\Delta }t)=(\mathrm{\Delta }t/2)𝒢_0(\mathrm{\Delta }t)(\mathrm{\Delta }t/2),$$
(5.11)
where $`𝒦(\mathrm{\Delta }t)`$ is regular in the vicinity of $`\mathrm{\Delta }t=0`$ and vanishes whenever $`_0`$ and $`𝒱`$ commute, and $`𝒢_0(\mathrm{\Delta }t)`$ is the amplification matrix when $`𝒱`$ is set to zero. The von Neumann condition is satisfied for $`𝒢_0(\mathrm{\Delta }t)`$ for a sufficiently small $`\mathrm{\Delta }t>0`$. Hence, powers of $`𝒢_0(\mathrm{\Delta }t)`$ applied to $`\mathrm{\Psi }(0)`$ cannot produce any exponential norm growth. Powers of $`𝒢_V(\mathrm{\Delta }t)`$ differ from those of $`𝒢_0(\mathrm{\Delta }t)`$ by factors that are powers of $`e^{\mathrm{\Delta }t𝒱}`$ and, hence, can only produce exponential attenuation of the amplitude. Indeed, let $`\mathrm{\Psi }_V(t)=e^{t𝒱}\mathrm{\Psi }(0)`$. Then $`/t(\mathrm{\Psi }_V,\mathrm{\Psi }_V)=2(\mathrm{\Psi }_V,𝒱\mathrm{\Psi }_V)0`$ since $`𝒱`$ is negative semidefinite. Thus, the approximate solution produced by the amplification matrix $`𝒢_V(\mathrm{\Delta }t)`$ has no exponential growth, while differing, relative to the energy norm, from that produced by $`𝒢_{}(\mathrm{\Delta }t)`$ by order of $`O(\mathrm{\Delta }t^2)`$. Therefore the modified leapfrog scheme can be made conditionally stable and as accurate as desired by reducing the time step.

It should be noted that our arguments do not prove that there cannot be any exponential growth of the norm $`\mathrm{\Psi }(t)=(\mathrm{\Psi }(t),\mathrm{\Psi }(t))^{1/2}`$ in the modified leapfrog scheme (5.9). All we can claim is that $`\mathrm{\Psi }(t)\mathrm{exp}(Kt)`$ which is also true for the conventional leapfrog scheme. The difference is that in the modified leapfrog scheme $`KO(\mathrm{\Delta }t^2)`$ as we have argued (a consequence of (5.11) and uniform boundedness of powers of $`𝒢_V`$), while in the conventional leapfrog scheme the constant $`K`$ is independent of $`\mathrm{\Delta }t`$. Hence a possible exponential growth cannot be suppressed by reducing the time step in (5.3), while it can be done in (5.9).

### V An example of the Lorentz model

To illustrate our general method we give an example of the Lorentz model commonly used to describe dispersive media. In the field representation of the Hamiltonian for the Lorentz model, we make the following decomposition

$$^F=\left(\begin{array}{cc}_0& 𝒱_{FM}\\ 𝒱_{MF}& 0\end{array}\right)+\left(\begin{array}{cc}0& 0\\ 0& _M^F\end{array}\right)_0^F+𝒱^F.$$
(5.12)
Substituting this decomposition into (5.9) we arrive at the following scheme

$`\psi ^F(t+\mathrm{\Delta }t)`$ $`=`$ $`\psi ^F(t\mathrm{\Delta }t)+2\mathrm{\Delta }t_0\psi ^F(t)+2\mathrm{\Delta }t{\displaystyle \underset{a}{}}𝒱_{FMa}\xi ^a(t),`$ (5.13)
$`\xi ^a(t+\mathrm{\Delta }t)`$ $`=`$ $`e^{2\mathrm{\Delta }t_{Ma}^F}\xi ^a(t\mathrm{\Delta }t)+2\mathrm{\Delta }te^{\mathrm{\Delta }t_{Ma}^F}𝒱_{MFa}\psi ^F(t),`$ (5.14)
$`e^{t_{Ma}^F}`$ $`=`$ $`e^{\gamma _at}\left[\mathrm{cosh}\stackrel{~}{\nu }_at+{\displaystyle \frac{\mathrm{sinh}\stackrel{~}{\nu }_at}{\stackrel{~}{\nu }_a}}\left(_{Ma}^F+\gamma _a\right)\right],`$ (5.15)
where $`\stackrel{~}{\nu }_a=(\gamma _a^2\omega _a^2)^{1/2}`$ and the six-dimensional columns $`\xi _a`$ are introduced in Section 3. The exponential (5.15) is easy to compute by expanding $`_{Ma}^F`$ in the Pauli matrix basis (a basis for the Lie algebra $`su(2)`$) and then by using the well known formula for the exponential of a linear combination of Pauli matrices. For small attenuation, $`\gamma _a<\omega _a`$, we get $`\stackrel{~}{\nu }_a=i\nu _a`$. The hyperbolic functions in (5.15) become the trigonometric ones and $`\stackrel{~}{\nu }_a`$ is replaced by $`\nu _a`$. Eigenvalues of the matter Hamiltonian are $`\lambda _a=\gamma _a\pm \stackrel{~}{\nu }_a`$. Hence, $`\mathrm{Re}\lambda _a<0`$ and amplitudes of the matter fields are always exponentially attenuated as $`t\mathrm{}`$, unless $`\gamma _a=0`$ leading to $`\mathrm{Re}\lambda _a=0`$.

The stability is ensured if $`_0^F`$ satisfies the von Neumann condition. Let $`k_{max}`$ be the maximal norm of all wave vectors of the initial wave packet and $`\omega _p^{max}`$ be the maximal value of $`\omega _p=(_a\omega _{pa}^2)^{1/2}`$ as a function of position, then a sufficient condition for stability reads

$$\mathrm{\Delta }t\sqrt{c^2k_{max}^2+(\omega _p^{max})^2}1.$$
(5.16)
Here the idea of the frozen coefficients  has been used again. The left hand side of inequality (5.16) is nothing but the spectral radius of $`\mathrm{\Delta }t_0^F`$ with frozen plasma frequencies so that $`\omega _p=\omega _p^{max}`$. Note that it is not difficult to solve the characteristic equation for $`_0^F`$ with frozen plasma frequencies by using the Fourier basis. The scheme (5.9) becomes especially simple in the case of small attenuation $`\gamma _a<\omega _a`$. In the complex representation of the auxiliary fields (3.8) (cf. (3.7)) the matter Hamiltonians $`_{Ma}^F`$ are diagonal and the action of its exponential is reduced to multiplication by a complex number $`e^{i\nu _a\mathrm{\Delta }t}`$.

Finally, it should be mentioned that, by rearranging operators in the split, namely, by moving $`𝒱_{FM}`$ and $`𝒱_{MF}`$ to $`𝒱^F`$ in (5.12), the stability condition (5.16) can be weakened to $`\mathrm{\Delta }tck_{max}1`$. This would come at the price of having a more complicated expression for $`(\mathrm{\Delta }t)`$. In the case of the Lorentz model it can still be computed analytically. The new split can also be viewed as the use of the induction representation in the modified leapfrog scheme, $`^I=_0^I+𝒱^I`$ where $`_0^I`$ contains only the blocks of $`^I`$ with the $`\mathbf{}`$ operator. The proof of the weaker stability condition can be found in .

## 6  Metal gratings in the Drude formalism

Here we apply our method to gratings made of a metal whose optical properties are described by the Drude formalism. This is an actual numerical scheme used in simulations in Section 7. The metal dielectric constant as a function of frequency is given by

$$\epsilon (\omega )=1+\stackrel{~}{\chi }(\omega )=1\frac{\omega _p^2}{\omega (\omega +i\eta )},$$

where $`\omega _p`$ is the plasma frequency, which is zero in the vacuum and constant in the region occupied by the metal as shown in Fig.1, and $`\eta `$ is the absorption. The model coincides with a one-resonant Lorentz model if $`\omega _0=0`$ and $`\eta 2\gamma `$. To satisfy the Gauss law exactly in simulations, we use the induction representation according to Section 5.II. An auxiliary field is chosen so that its first component equals $`𝐏`$ and the second is denoted $`𝝃`$. The Hamiltonian evolution equations are taken in the form

$`\dot{𝐃}`$ $`=`$ $`c\mathbf{}\times 𝐁,`$
$`\dot{𝐁}`$ $`=`$ $`c\mathbf{}\times \left(𝐃𝐏\right),`$
$`\dot{𝐏}`$ $`=`$ $`\eta 𝝃,`$
$`\dot{𝝃}`$ $`=`$ $`\eta 𝝃{\displaystyle \frac{\omega _p^2}{\eta }}\left(𝐃𝐏\right).`$
To apply the modified leapfrog scheme the Hamiltonian is split into the sum

$`^I`$ $`=`$ $`\left(\begin{array}{cccc}0& c\mathbf{}\times & 0& 0\\ c\mathbf{}\times & 0& c\mathbf{}\times & 0\\ 0& 0& 0& \eta \\ \omega _p^2\eta ^1& 0& \omega _p^2\eta ^1& \eta \end{array}\right)=_0^I+𝒱^I,`$ (6.1)
$`𝒱^I`$ $`=`$ $`\mathrm{𝖽𝗂𝖺𝗀}(0,0,0,\eta ).`$
Clearly, $`𝒱^I`$ is negative semidefinite because $`\eta >0`$. The stability of the modified leapfrog scheme requires that eigenvalues of $`_0^I`$ have zero real parts. This is indeed the case. In the Fourier basis, the $`\mathbf{}`$ operator becomes $`i𝐤`$. It is straightforward to find the characteristic polynomial $`det(_0^I\lambda )`$. Its nonzero roots are $`\lambda =\pm i\sqrt{c^2𝐤^2+\omega _p^2}`$. Hence the scheme is stable if the time step is chosen so that the condition (5.16) is satisfied. In particular, $`k_{\mathrm{max}}`$ can be set to $`k_0`$ being the step of the reciprocal lattice and $`\omega _p^{max}`$ is the plasma frequency of the metal (silver in our simulations, see Section 7). An explicit form of the modified leapfrog scheme for the split (6.1) reads

$`𝐃(t+\mathrm{\Delta }t)`$ $`=`$ $`𝐃(t\mathrm{\Delta }t)+2c\mathrm{\Delta }t\mathbf{}\times 𝐁(t),`$
$`𝐁(t+\mathrm{\Delta }t)`$ $`=`$ $`𝐁(t\mathrm{\Delta }t)2c\mathrm{\Delta }t\mathbf{}\times \left[𝐃(t)𝐏(t)\right],`$
$`𝐏(t+\mathrm{\Delta }t)`$ $`=`$ $`𝐏(t\mathrm{\Delta }t)2\eta \mathrm{\Delta }t𝝃(t),`$
$`𝝃(t+\mathrm{\Delta }t)`$ $`=`$ $`e^{2\eta \mathrm{\Delta }t}𝝃(t\mathrm{\Delta }t)2\mathrm{\Delta }t\omega _p^2\eta ^1e^{\eta \mathrm{\Delta }t}\left[𝐃(t)𝐏(t)\right].`$
The action of the curl is computed by the fast Fourier method in combination with a change of variables that enhances the sampling efficiency near the metal-vacuum interface. The details are in Section 7.

## 7 Results for extraordinary transmission gratings

To test our method we applied it to study transmission properties of metal and dielectric gratings suspended in vacuum. Transmission and reflection gratings have been subject of a number of experimental and theoretical works , . The interest is stimulated by nearly 100% transmission or reflection within narrow wavelength range with possibility of using such grating as efficient filters. Moreover, extraordinary optical transmission has been observed in the $`2D`$ hole arrays , further stimulating theoretical and experimental interest to transmission properties of nanostructured materials .

In Fig.1 we schematically represent our system comprising a metallic or dielectric slab with gratings. $`D`$ is a grating period, $`a`$ the size of the grating and $`h`$ its thickness. For the sake of comparison with previously published results we have chosen $`D=1.75\mu m`$ and $`a=0.30\mu m`$. Transmission and reflection coefficients are computed as a function of the grating thickness $`h`$. An incident electromagnetic wave packet propagates along the $`z`$-direction, normal to the slab. The polarization is such that the electric field vector is perpendicular to the gratings, while the magnetic field is parallel to them (the so called $`p`$ -polarization). Calculations are performed in a finite $`(x,z)`$ box of the size -$`15D<z<12D`$, $`D/2<x<D/2`$. A uniform mesh of typically $`256`$ knots is used in the $`x`$-direction and a nonuniform mesh of $`512`$ knots generated by the change of variables is used in the $`z`$-direction. The change of variables is used to enhance the sampling efficiency near the two interfaces in the $`z`$-direction. Periodic boundary conditions are insured in $`x`$ through the pseudospectral approach based on the Fast Fourier Transform. To suppress an artificial reflection of the wave packet, absorbing layers are introduced at the box boundaries $`z=\pm 15D`$. The initial wave packet is Gaussian, which allows us to obtain transmission (reflection) coefficients within the entire frequency bandwidth of the initial wave packet by a single simulation. In what follows we are mainly interested in transmission (reflection) of the radiation with wavelengths larger than the grating period (zero order diffraction). Thus, reflected or transmitted waves propagate in the direction of the $`z`$-axis, the same as the incident radiation. Note, however, that this is not a limitation of our method which allows for a priori extraction of the entire scattering matrix for all wave vectors.

In Fig.2 we show transmission and reflection coefficients obtained for metallic gratings of variable depth. The dielectric response of the metal is described within the Drude formalism (see Section 6). Following  we use $`\omega _p=9eV`$ and $`\eta =0.1eV`$ representative for silver. As clearly seen in the Figure, the transmission coefficients exhibit narrow resonances for certain wavelengths. With increasing thickness of the gratings, $`h`$, the number of resonant structures increases. Our results are in a full agreement with previously published theoretical studies . The only difference being, that the model used in  assumes perfect metal surfaces inside the gratings, neglecting possible absorption of the radiation. This leads to $`100\%`$ transmission at the resonant frequencies. In our case, a part of radiation is absorbed by the metal so that the transmission never reaches $`100\%`$.

The observed structures in transmission are associated with resonant modes of the electromagnetic field produced by coupled surface electromagnetic modes (called surface plasmon polaritons) and waveguide modes inside the gratings. Some of these resonances posses relatively long lifetimes. This can be immediately inferred from their width, as, e.g., the resonance located at $`\lambda =1.1D`$ in the case of the grating with thickness $`h=1.4\mu m`$. Another way to observe the trapped field resonances is to look at the time dependence of the field transmitted in the $`z`$-direction. The signal in Fig. 3 is registered by a detector placed at the distance of $`3.5D`$ behind the gratings. As seen in the figure, as soon as the resonances are populated by an incident $`25fs`$-long pulse, they radiate the field during at least $`125fs`$. Here the radiation time is determined by the lifetime of the narrowest resonance located at $`\lambda =1.1D`$. It is worth mentioning that, as shows the sum of the transmission and reflection coefficients, the total absorption is largest at the resonance positions, i.e., when the interaction time between the radiation (trapped mode) and the metal is large.

Finally, Fig. 4 shows the field structure of a trapped mode corresponding to the narrowest resonance observed with the $`h=2.4\mu m`$ thickness grating. The $`E_x`$ component of the electric field is presented. It is obtained via sufficiently long time propagation so that contributions from less long-lived states vanish.

In Figs. 5 and 6 we show reflection coefficients for dielectric gratings suspended in vacuum. Simulations are performed using the conventional leapfrog scheme (5.3) applied to (4.5). Since the attenuation is absent, the scheme is stable if the time step satisfies (5.6). The dielectric material is modeled through the frequency independent dielectric constant $`\epsilon =2`$ (Fig.5) and $`\epsilon =4`$ (Fig.6). Only the reflection coefficient is shown here since there is no absorption of the radiation in dielectric so that the transmission can simply be inferred from the unitarity of the scattering matrix. Without gratings the dielectric slabs are basically transparent in both cases. Introducing grating structures results in a complete reflection of the incident radiation within an extremely narrow wavelength bandwidth. The associated guided mode resonances have been extensively discussed in the literature . In the case of $`\epsilon =2`$ the resonances are so narrow that extraction of frequency dependent transmission and reflection coefficients by the Fourier transform of the scattered wave would require too large propagation time (see also Fig. 7). We had to stop our wave packet propagation before the radiation emitted by resonances ended. This explains why the reflection coefficients do not reach its maximal value $`1`$ in this case. Consistently with the metal case, the number of resonant structures increases with increase of the width $`h`$ of the gratings. The width of resonances increases with $`\epsilon `$ as follows from the comparison of Figs. 5 and 6. Note also that resonances are associated with Fano profiles that usually arise because of the interference between the non-resonant and resonant contributions to the scattered wave. Such narrow reflection structures in the case of dielectric gratings have been usually studied by stationary methods in the frequency domain. An important advantage of the time-dependent study is that one has an immediate access to all the details the temporal evolution of electromagnetic fields in any desired part of the system.

In Fig. 7 we show the time dependence of the field transmitted in the $`z`$-direction for gratings characterized by $`\epsilon =2`$ and $`h=0.8\mu m`$. With these parameters there is only one resonance in the reflection spectrum. The signal is registered by a detector placed at the distance of $`3.5D`$ behind the gratings. First, we observe that a $`25fs`$ incident pulse is transmitted through the structure without modification. The lasing effect, when the transmitted field is followed by basically monochromatic radiation, can clearly be seen. For readability reasons we could not show the complete time evolution in the figure, but the lasing effect lasts for at least (!) $`2ps`$ reflecting an extraordinary long lifetime of the trapped resonant field. It is this radiation which comes in the phase opposite to the corresponding harmonic in the transmitted initial pulse and leads, finally, to the zero transmission at the corresponding frequency. The same lasing effect to the left from the grating structure is responsible for the $`100\%`$ reflection at the same frequency.

In Fig. 8 we show a typical structure of the field corresponding to a trapped (resonant) mode. The $`D_x`$ component of the electric induction is represented in the case of dielectric gratings with $`\epsilon =4.`$ The thickness of the dielectric slab is $`h=0.6\mu m`$. Note that in contrast to the metal grating structure, the field in the present case occupies the entire slab and not only the vacuum part of the grating.

## 8 Conclusions

We have developed a time domain algorithm for the initial value problem for the Maxwell’s theory of linear passive media. The algorithm is based on (i) the Hamiltonian formalism for evolution differential equations, (ii) the Fourier pseudospectral method in which the sampling efficiency in designated space regions is enhanced by a suitable change of variables, and (iii) the modified leapfrog scheme. We have analyzed the stability of the algorithm and found explicit stability conditions when passive media are described by multi-resonant Lorentz models. We have implemented and tested our algorithm for extraordinary transmission and reflection gratings whose optical properties have been studied in a number of theoretical and experimental works. Numerical simulations based on our algorithm are shown to produce extremely accurate data for the well studied far-field (zero-order diffraction) at relatively low computational costs. A single simulation of an incident wideband wave packet is sufficient to determine transmission and reflection properties of the gratings in the frequency range of the initial wave packet. In addition, our algorithm allows us to see a real time dynamics of formation of long-living resonant excitations of electromagnetic fields in the grating as well as formation of transmitted and/or reflected wave fronts in the entire frequency range of the initial wave packet. It is believed that our algorithm would be useful for numerical studies of other nanostructured materials.

Acknowledgments

S.V.S thanks the Laboratory of Atomic and Molecular Collisions of the University of Paris-Sud for the warm hospitality and Dr. R. Albanese (Brooks, Air Force Base, TX) for his continued interest and support of this project. The work of S.V.S. has been supported in part by US Air Force Grants F4920-03-10414 and F49620-01-1-0473.

Figures

Fig. 1 A schematic representation of the studied system. The incident wave packet propagates along the normal to the slab containing the gratings (z-direction).

Fig. 2. Calculated zero-order transmission (solid lines) and reflection (dashed lines) coefficients for metallic gratings described in the text. The results are presented as a function of the wavelength $`\lambda `$ of the incident radiation measured in the units of the period of the gratings, $`D`$. Different panels of the figure correspond to the different thickness $`h`$ of the gratings, as indicated.

Fig. 3. The electric field measured by a detector placed behind the metallic gratings with thickness $`h=1.4\mu m`$. Only the field corresponding to the zero-order transmitted wave propagating along the $`z`$-axis is represented. It is obtained by the Fourier analysis of the $`x`$-coordinate dependence of the field at the detector position. The signal is shown as a function of time measured in femtoseconds.

Fig. 4. A snapshot of the $`x`$-component of the electric field, $`E_x`$ for the grating thickness $`h=2.4\mu m`$. The results are presented as a function of $`x`$ and $`z`$ coordinates measured in units of $`D`$. Red and blue colors correspond, respectively, to positive and negative values of the field. The snapshot has been produced after a sufficiently long propagation time so that the field pictured in it indeed corresponds to the long-lived resonance giving enhanced transmission at $`\lambda =1.15D`$ (see Fig.2).

Fig. 5. Calculated zero-order reflection coefficient for dielectric gratings with $`\epsilon =2`$. The results are presented as a function of the wavelength $`\lambda `$ of the incident radiation measured in the units of the period of the gratings, $`D`$. Different panels of the figure correspond to the different thickness $`h`$ of the gratings, as indicated.

Fig. 6. Calculated zero-order reflection coefficient for dielectric gratings with $`\epsilon =4`$. Results are presented as a function of the wavelength $`\lambda `$ of the incident radiation measured in the units of the period of the gratings, $`D`$. Different panels of the figure correspond to the different thickness $`h`$ of the gratings, as indicated.

Fig. 7. The electric field measured by a detector placed behind the dielectric grating with thickness $`h=0.8\mu m`$ and dielectric constant $`\epsilon =4`$. Only the field corresponding to the zero-order transmitted wave propagating along the $`z`$-axis is represented. It is obtained by the Fourier analysis of the $`x`$-coordinate dependence of the field at the detector position. The signal is shown as a function of time measured in femtoseconds.

Fig. 8. A snapshot of the $`x`$-component of the electric induction, $`D_x`$ for the dielectric grating with dielectric constant $`\epsilon =4`$ and thickness $`h=0.6\mu m`$. The results are presented as a function of $`x`$ and $`z`$ coordinates measured in units of $`D`$. Red and blue colors correspond, respectively, to positive and negative values of the induction. The snapshot has been produced after a sufficiently long propagation time so that the induction pictured in it indeed corresponds to the resonance giving enhanced transmission at $`\lambda 1.3D`$ (see Fig.6).