Alternative expression of the Bloch wave group velocity in loss-less periodic media using the electromagnetic field energy

Abstract

In periodic optical media, the group velocity is defined as the gradient with respect to wave-vector of the corresponding Bloch mode frequency dispersion curve, forming the photonic band structure. Instead of deducing it from the numerically computed photonic crystal band structure, the group velocity can be calculated directly from the integral of the Poynting vector over the crystal unit cell, the physical meaning of which is immediately perceivable. The related formula, which can be regarded as the application of Hellmann–Feynman theorem to electromagnetism, has been reported previously though without proof. We provide hereafter a full derivation of that formula starting from Maxwell’s equations and we discuss its usefulness in photonics.

Keywords:

• Photonic crystals
• group velocity
• slow light

1. Introduction

In a perfect crystal, electrons can travel long distances while being scattered by the periodic potential of the atoms, a property that is easily explained using the concept of Bloch waves [1]. In the context of optics, Bloch waves originate from coherent light scattering by regularly ordered and refractive-index-contrasted dielectric regions that form a so-called photonic crystal, a terminology coined in the early 90s [2]. Bloch modes are eigen mode solutions of Maxwell’s equations in the periodic potential defined by the spatially varying permittivity $\mathit{\epsilon }(\mathbf{r})$. They are labelled according to the wave-vector $\mathbf{k}$ and the band index n. Their dispersion relations, ${\mathit{\omega }}_n(\mathbf{k},n)$ ($\mathit{\omega }$: angular frequency), form the photonic band structure of the infinite periodic medium [2]. In three-dimensional lattices, it is not obvious to figure out intuitively how waves actually propagate in an arbitrary direction. The Poynting vector [3], when it is calculated at any point of the crystal unit cell, allows us to describe the flow of electromagnetic energy and to draw analogies with fluid mechanics (cf. Poynting’s theorem). Poynting vector is a local quantity, the spatial distribution of which needs to be calculated only over the unit cell due to periodic boundary conditions.

In a periodic medium, the group velocity for a given Bloch wave, defined as the gradient with respect to $\mathbf{k}$ of the mode frequency ${\mathit{\omega }}_n$: ${\mathbf{v}}_{\mathbf{g}}\equiv {\mathbf{\nabla }}_{\mathbf{k}}{\mathit{\omega }}_n$, is a global quantity, i.e. a constant vector. It does not depend on the coordinates but depends on the wave-vector $\mathbf{k}$ and the band index n. For a given mode propagation direction, it can be determined from the band structure by numerically computing the slope of the dispersion curve ${\mathit{\omega }}_n(\mathbf{k})$ associated with that mode, at the corresponding wave-vector position. An alternative to calculate the group velocity is to use ${\mathbf{v}}_g=\frac{1}{w}\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})$, where $\mathbf{S}(\mathbf{r})$ is the Poynting vector at arbitrary point $\mathbf{r}$ in the unit cell (the integration being performed over the unit cell volume) and w is the electromagnetic energy density per unit cell (Equation (13(13) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})}{{\mathcal{U}}_H+{\mathcal{U}}_E}\text{.}\end{array}$$(13) ) of this article). This formula, which expresses the equivalence between group velocity and energy velocity in a loss-less medium [4,5], tells us that the group velocity for a given mode is proportional to the integral of the Poynting vector of that mode over the unit cell. It can, therefore, be calculated in a straightforward manner from a local quantity, the Poynting vector, the physical meaning of which is immediately perceivable. In order to apply the above formula to a given periodic structure, the electric and magnetic fields must be known first at any point in the unit cell. In practice, $\mathbf{E}(\mathbf{r})$ and $\mathbf{H}(\mathbf{r})$ can be computed using numerical methods that solve Maxwell’s equations in periodic media (see for instance [6]). Then, knowing the fields everywhere, the integral of the Poynting vector and the energy density per unit cell can be calculated directly, leading to the group velocity. Although the formula, with some explanations, can be found in probably the most cited textbook on the theory of photonic crystals [7], no formal proof has been given to the best of our knowledge. It is the purpose of our article to give a proof of this formula.

The alternative expression of the Bloch wave group velocity is expected to be useful in research topics where the control of light wave group velocity matters, such as in applications related to so-called slow light [8,9]. In photocatalysis, for instance, inverse-opal-type 3D photonic crystals are developed [10] in order to enhance solar light absorption, hence photocatalytic activity, through the increase of the photon lifetime (decrease of group velocity) at frequencies located at the edges of the photonic band gap [11]. Relying on the above formula, photonic crystals can be designed in order to optimize concomitantly the flow of light (Poynting vector field pattern) and the group velocity of the slow photons in the inverse opal structure. Since the objective is light harvesting enhancement in the photoactive (i.e. optically absorbing) material, the medium permittivity can no longer be considered as real quantity, which precludes in principle the exploitation of the results of the present theory. However, the slow-photon-induced enhancement is targeted in spectral regions where material absorption is small [11] so that, to first approximation, the imaginary part of the permittivity can be considered as a small perturbation. Extending the results of this article to the case of a medium with complex permittivity is an interesting perspective to explore in the future.

The article is organized as follows. First, the eigenvalue problem formulation of light wave propagation in optical periodic media is recalled. Then, the expression of the Poynting vector is derived for Bloch waves and the group velocity formula is presented with a brief comment on the premises to its derivation. Next, the linear differential tensor operator acting on Bloch modes is defined and its expression is derived. Finally, a proof of the formula is provided. The last section gives a conclusion. Details about mathematical formalism and calculation rules used in this article are given in Appendix 1.

2. Eigenvalue problem formulation of light wave propagation in periodic media

Starting from Maxwell’s equations, the propagation of light wave in an arbitrary medium can be reformulated as an eigenvalue problem [7]. This is the approach we follow here to describe propagation in a periodic medium. In the more general case of an inhomogeneous medium, which is assumed to be non-magnetic, loss-less and isotropic, the wave equation is written as:(1) $$\begin{array}{c}\hfill \mathrm{\nabla }\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}\mathrm{\nabla }\times \mathbf{H}(\mathbf{r})\right)={\left(\frac{\mathit{\omega }}{c}\right)}^2\mathbf{H}(\mathbf{r})\text{,}\end{array}$$(1)

where $\mathbf{H}(\mathbf{r})$ is the magnetic field (time dependence $e^{-i\mathit{\omega }t}$) and $\mathit{\epsilon }(\mathbf{r})$ is the medium permittivity (real and scalar quantity). This master equation is subject to the transversality condition: $\mathrm{\nabla }·\mathbf{H}(\mathbf{r})=0$.

The problem of solving (1(1) $$\begin{array}{c}\hfill \mathrm{\nabla }\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}\mathrm{\nabla }\times \mathbf{H}(\mathbf{r})\right)={\left(\frac{\mathit{\omega }}{c}\right)}^2\mathbf{H}(\mathbf{r})\text{,}\end{array}$$(1) ) for the harmonic modes $\mathbf{H}(\mathbf{r})$ can be recast into an eigenvalue problem by defining the linear differential vector operator [7](2) $$\begin{array}{c}\hfill \widehat{\mathrm{\Theta }}\equiv \mathrm{\nabla }\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}\mathrm{\nabla }\times \right)\text{,}\end{array}$$(2)

so that(3) $$\begin{array}{c}\hfill \widehat{\mathrm{\Theta }}\mathbf{H}(\mathbf{r})={\left(\frac{\mathit{\omega }}{c}\right)}^2\mathbf{H}(\mathbf{r})\text{,}\end{array}$$(3)

where the eigenvectors $\mathbf{H}(\mathbf{r})$ are the spatial patterns of the harmonic modes and the eigenvalues ${\left(\frac{\mathit{\omega }}{c}\right)}^2$ are proportional to the squared mode frequencies. It can be shown that the operator $\widehat{\mathrm{\Theta }}$ is Hermitian, i.e. the mode frequency to the second power $(\mathbf{F},\widehat{\mathrm{\Theta }}\mathbf{G})=(\widehat{\mathrm{\Theta }}\mathbf{F},\mathbf{G})$ where the parentheses denote the inner product [7].

Now, if the inhomogeneous medium is assumed to be a periodic arrangement of two loss-less dielectric constituents (${\mathit{\epsilon }}_1,{\mathit{\epsilon }}_2$ real) in the three dimensions of space, i.e. 3D photonic crystal, the modes take the form of Bloch states [7](4) $$\begin{array}{c}\hfill {\mathbf{H}}_{\mathbf{k}}(\mathbf{r})={\mathbf{u}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}\text{,}\end{array}$$(4)

where $\mathbf{k}={\sum }_j{k}_j{\mathbf{b}}_j$ (${\mathbf{b}}_j$: primitive reciprocal lattice vectors) is the Bloch wave-vector and ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})={\mathbf{u}}_{\mathbf{k}}(\mathbf{r}+\mathbf{R})$ is a periodic function (vector field) for all lattice vectors $\mathbf{R}=l{\mathbf{a}}_1+m{\mathbf{a}}_2+n{\mathbf{a}}_3$ (${\mathbf{a}}_j$: primitive lattice vectors; l, m, n: integers). Similarly, the Bloch states for the electric field are noted(5) $$\begin{array}{c}\hfill {\mathbf{E}}_{\mathbf{k}}(\mathbf{r})={\mathbf{e}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}\text{.}\end{array}$$(5)

Note that $\mathrm{\nabla }·{\mathbf{H}}_{\mathbf{k}}(\mathbf{r})=0$, i.e. the transversality condition, leads to $(i\mathbf{k}+\mathrm{\nabla })·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})=0$, hence(6) $$\begin{array}{c}\hfill \mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})=-i\mathbf{k}·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{.}\end{array}$$(6)

From Maxwell’s equation $\mathrm{\nabla }\times {\mathbf{H}}_{\mathbf{k}}(\mathbf{r})=-i\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r}){\mathbf{E}}_{\mathbf{k}}(\mathbf{r})$, we also deduce(7) $$\begin{array}{c}\hfill {\mathbf{e}}_{\mathbf{k}}(\mathbf{r})=\frac{i}{\mathit{\omega }{\mathit{\epsilon }}_0}\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{.}\end{array}$$(7)

Inserting the Bloch state (4(4) $$\begin{array}{c}\hfill {\mathbf{H}}_{\mathbf{k}}(\mathbf{r})={\mathbf{u}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}\text{,}\end{array}$$(4) ) into the master Equation (3(3) $$\begin{array}{c}\hfill \widehat{\mathrm{\Theta }}\mathbf{H}(\mathbf{r})={\left(\frac{\mathit{\omega }}{c}\right)}^2\mathbf{H}(\mathbf{r})\text{,}\end{array}$$(3) ) and applying twice $\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}=e^{i\mathbf{k}·\mathbf{r}}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k$ leads to another eigenvalue equation(8) $$\begin{array}{c}\hfill {\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})={\left(\frac{\mathit{\omega }(\mathbf{k})}{c}\right)}^2{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{,}\end{array}$$(8)

where ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ is a new Hermitian operator defined by(9) $$\begin{array}{c}\hfill {\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\equiv (i\mathbf{k}+\mathrm{\nabla })\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times \right)\text{.}\end{array}$$(9)

The solutions of this eigenvalue equation take the form of an infinite discrete set of modes with frequencies ${\mathit{\omega }}_n(\mathbf{k})$, i.e. the photonic band structure of the photonic crystal [7].

3. Poynting vector for Bloch waves

Let us now derive the general expression of the Poynting vector for Bloch waves in periodic media. The Poynting vector $\mathbf{S}$ at an arbitrary point of space $\mathbf{r}$ can be written as:(10) $$\begin{array}{c}\hfill \mathbf{S}(\mathbf{r})=\frac{1}{2}\mathfrak{R}\left[{\mathbf{E}}_{\mathbf{k}}(\mathbf{r})\times {\mathbf{H}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]=\frac{1}{2}\mathfrak{R}\left[{\mathbf{e}}_{\mathbf{k}}(\mathbf{r})\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\text{,}\end{array}$$(10)

where Equations (4(4) $$\begin{array}{c}\hfill {\mathbf{H}}_{\mathbf{k}}(\mathbf{r})={\mathbf{u}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}\text{,}\end{array}$$(4) )–(5(5) $$\begin{array}{c}\hfill {\mathbf{E}}_{\mathbf{k}}(\mathbf{r})={\mathbf{e}}_{\mathbf{k}}(\mathbf{r})e^{i\mathbf{k}·\mathbf{r}}\text{.}\end{array}$$(5) ) were used to obtain the second equality. By substituting (7(7) $$\begin{array}{c}\hfill {\mathbf{e}}_{\mathbf{k}}(\mathbf{r})=\frac{i}{\mathit{\omega }{\mathit{\epsilon }}_0}\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{.}\end{array}$$(7) ) for ${\mathbf{e}}_k$, we find$$\begin{array}{cc}\hfill \mathbf{S}(\mathbf{r})& =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\mathfrak{R}\left[i(i\mathbf{k}\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right.\hfill \\ \hfill & \left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\hfill \\ \hfill & =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\mathfrak{R}\left[i({\mathbf{u}}_{\mathbf{k}}(\mathbf{r})(i\mathbf{k}·{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r}))-i\mathbf{k}|{\mathbf{u}}_{\mathbf{k}}(\mathbf{r}){|}^2)\right.\hfill \\ \hfill & \left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\hfill \\ \hfill & =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\mathfrak{R}\left[i{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})(\mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r}))+\mathbf{k}{|{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})|}^2\right.\hfill \\ \hfill & \left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\text{,}\hfill \end{array}$$

where we used $(\mathbf{a}\times \mathbf{b})\times \mathbf{c}=\mathbf{b}(\mathbf{a}·\mathbf{c})-\mathbf{a}(\mathbf{b}·\mathbf{c})$ and the conjugate of (6(6) $$\begin{array}{c}\hfill \mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})=-i\mathbf{k}·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{.}\end{array}$$(6) ) in the second and third equalities, respectively. Since $\mathbf{k}$ is a real quantity, we can rewrite this result as,(11) $$\begin{array}{cc}\hfill \mathbf{S}(\mathbf{r})& =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\left(\mathbf{k}|{\mathbf{u}}_{\mathbf{k}}{(\mathbf{r})|}^2\right.\hfill \\ \hfill & +\mathfrak{R}\left[i{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})(\mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r}))\right.\hfill \\ \hfill & \left.\left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\right)\text{.}\hfill \end{array}$$(11)

In the case of a homogeneous medium ($\mathit{\epsilon }(\mathbf{r})\equiv \mathit{\epsilon }$) and for a plane wave (${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\equiv {\mathbf{H}}_0$), the first and the second terms inside the square brackets vanish and we find the well-known result that the Poynting vector is uniform in space and aligned along the wave propagation direction $\mathbf{k}$, i.e. $\mathbf{S}=\frac{|{\mathbf{H}}_0{|}^2}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }}\mathbf{k}$. This is no longer the case for a Bloch wave in a periodic medium, for which both the norm and the orientation of the Poynting vector vary in space and are related to the local value of ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})$. We note that both the first and the second terms inside the square brackets take complex values, only the real parts of them being kept: $\mathfrak{R}\left[i{\mathbf{u}}_k(\mathrm{\nabla }·{\mathbf{u}}_k^{\ast })\right]=\mathfrak{R}{\mathbf{u}}_k(\mathrm{\nabla }·\mathfrak{I}{\mathbf{u}}_k)-\mathfrak{I}{\mathbf{u}}_k(\mathrm{\nabla }·\mathfrak{R}{\mathbf{u}}_k)$ and $\mathfrak{R}\left[i(\mathrm{\nabla }\times {\mathbf{u}}_k^{\ast })\times {\mathbf{u}}_k\right]=(\mathrm{\nabla }\times \mathfrak{R}{\mathbf{u}}_k)\times \mathfrak{I}{\mathbf{u}}_k-(\mathrm{\nabla }\times \mathfrak{I}{\mathbf{u}}_k)\times \mathfrak{R}{\mathbf{u}}_k$.

Even though ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r})$ and related quantities are subjected to periodic boundary conditions, no simplification appears when integrating the Poynting vector over the volume of the photonic crystal unit cell.

However, in the particular case of a one-dimensional (1D) photonic crystal, so called Bragg layer stack, it is noteworthy that $\mathbf{S}$ and $\mathbf{k}$ are colinear and a simplified expression can be obtained for the spatial evolution of the projection of $\mathbf{S}$ on the direction (z-axis) perpendicular to the planar interfaces, i.e. $S_z(z)$ [12].

4. Group velocity formula

Contrary to the case of a homogenous isotropic medium (and also 1D photonic crystal) where the electromagnetic energy carried by the plane wave flows along the wave-vector direction , the direction and the speed of the Bloch wave in a 2D or 3D periodic medium are given by the group velocity, which depends on both the wave-vector $\mathbf{k}$ and the photonic band number n: ${\mathbf{v}}_g\equiv {\mathbf{\nabla }}_{\mathbf{k}}{\mathit{\omega }}_n$, where ${\mathbf{\nabla }}_{\mathbf{k}}$ is the gradient with respect to $\mathbf{k}$.

Using the formalism developed above, an expression of the group velocity in terms of eigenvectors ${\mathbf{u}}_k$ can be derived [7]. We first differentiate the eigenvalue Equation (8(8) $$\begin{array}{c}\hfill {\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})={\left(\frac{\mathit{\omega }(\mathbf{k})}{c}\right)}^2{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{,}\end{array}$$(8) ) for a given band index n with respect to $\mathbf{k}$ and then take the inner product (see Notes) with ${\mathbf{u}}_k$ on both sides:$$\left({\mathbf{u}}_k,{\mathbf{\nabla }}_{\mathbf{k}}\left[{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_{\mathbf{k}}\right]\right)=\left({\mathbf{u}}_k,{\mathbf{\nabla }}_{\mathbf{k}}\left[{\left(\frac{\mathit{\omega }}{c}\right)}^2{\mathbf{u}}_k\right]\right)\text{.}$$

By differentiating the products ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_{\mathbf{k}}$ and ${\left(\frac{\mathit{\omega }}{c}\right)}^2{\mathbf{u}}_k$, we obtain$$\begin{array}{cc}& \left({\mathbf{u}}_k,\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k+{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{u}}_k\right)\hfill \\ \hfill & =\left({\mathbf{u}}_k,2\frac{\mathit{\omega }}{c^2}{\mathbf{v}}_g{\mathbf{u}}_k+{\left(\frac{\mathit{\omega }}{c}\right)}^2{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{u}}_k\right)\text{.}\hfill \end{array}$$

Since ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ is Hermitian, we have $\left({\mathbf{u}}_k,{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{u}}_k\right)=\left({\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k,{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{u}}_k\right)$ and the ${\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{u}}_k$ terms cancel out on both sides after using (8(8) $$\begin{array}{c}\hfill {\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})={\left(\frac{\mathit{\omega }(\mathbf{k})}{c}\right)}^2{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{,}\end{array}$$(8) ). The remaining terms lead to the following expression of the group velocity [7]:(12) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{c^2}{2\mathit{\omega }}\frac{\left({\mathbf{u}}_k,\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k\right)}{\left({\mathbf{u}}_k,{\mathbf{u}}_k\right)}\text{.}\end{array}$$(12)

This formula can be regarded as the result of the application of Hellmann–Feynman theorem [13] to electromagnetism. This theorem relates variations in an operator involved in an inner product to corresponding variations in the operator eigenvalue. It shows that, in order to compute the derivative of an eigenvalue (here ${\mathit{\omega }}_n$) with respect to a parameter of the operator (here $\mathbf{k}$), one needs only to know an eigenvector (here ${\mathbf{u}}_k$) and the derivative of the operator. The tensor operator $\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]$ is defined from the differentiation of the vector operator ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ with respect to $\mathbf{k}$. An expression of this tensor, which is not given explicitly in [7], will be derived in the next section.

The denominator is a scalar quantity equal to $\left({\mathbf{u}}_k,{\mathbf{u}}_k\right)=\int d^3\mathbf{r}|{\mathbf{u}}_{\mathbf{k}}{(\mathbf{r})|}^2=\int d^3\mathbf{r}{|{\mathbf{H}}_{\mathbf{k}}(\mathbf{r})|}^2=\frac{4}{{\mathit{\mu }}_0}{\mathcal{U}}_H=\frac{2}{{\mathit{\mu }}_0}({\mathcal{U}}_H+{\mathcal{U}}_E)$, i.e. it is proportional to the sum of the electric and magnetic energy contained in a unit cell.

The numerator must, therefore, be a real vector in order to be consistent with the definition of ${\mathbf{v}}_g$. We will demonstrate hereafter that$$\frac{c^2}{2\mathit{\omega }}\left({\mathbf{u}}_k,\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k\right)=\frac{2}{{\mathit{\mu }}_0}\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})\text{,}$$

so that, all together, we have(13) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})}{{\mathcal{U}}_H+{\mathcal{U}}_E}\text{.}\end{array}$$(13)

The above quantity, where the Poynting vector at an arbitrary point of the crystal unit cell is given by (11(11) $$\begin{array}{cc}\hfill \mathbf{S}(\mathbf{r})& =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\left(\mathbf{k}|{\mathbf{u}}_{\mathbf{k}}{(\mathbf{r})|}^2\right.\hfill \\ \hfill & +\mathfrak{R}\left[i{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})(\mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r}))\right.\hfill \\ \hfill & \left.\left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\right)\text{.}\hfill \end{array}$$(11) ), is the velocity of the energy propagation in the photonic crystal medium for a given mode. More precisely stated, in a loss-less (i.e. non-absorbing) periodic medium and in the absence of sources, we can write a continuity equation $⟨\mathbf{S}⟩=⟨w⟩{\mathbf{v}}_g$ where $⟨\cdots ⟩$ denotes the average of a spatial quantity over a unit cell and w is the electromagnetic energy density. Equation (13(13) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})}{{\mathcal{U}}_H+{\mathcal{U}}_E}\text{.}\end{array}$$(13) ) is the central result of this article, a rigorous demonstration of which follows.

5. Tensor operator acting on Bloch modes

We define the tensorial operator ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ by analogy with the definition of the Jacobian matrix of a vector field:(14) $$\begin{array}{c}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]=\left(\begin{array}{ccc}\frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}\\ \frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}\\ \frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}\end{array}\right).\end{array}$$(14)

Note that the components of the gradient of each component ${\mathbf{e}}_i·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ of the vector operator ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ appear along the columns of the matrix representing the tensor operator. In a shorthand notation where comma separates columns of the matrix: $[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}]=({\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}},{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}},{\mathbf{\nabla }}_{\mathbf{k}}{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}})$.

Application of this tensorial operator to the Bloch mode produces a second-rank tensor which can be written as:(15) $$\begin{array}{c}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k=\sum _{i,j}\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_j·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\right){\mathbf{e}}_i{\mathbf{e}}_j\text{,}\end{array}$$(15)

where the dyadic notation ${\mathbf{e}}_i{\mathbf{e}}_j$ is used and i, j are row and column indices in the above matrix representation.

A physical interpretation can be given for this second-rank tensor. Indeed, if we multiply the transpose of the matrix (15(15) $$\begin{array}{c}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k=\sum _{i,j}\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_j·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\right){\mathbf{e}}_i{\mathbf{e}}_j\text{,}\end{array}$$(15) ) by a unit vector $\mathbf{n}$ pointing in an arbitrary direction, we obtain a vector the components of which are the projections, in that direction, of the gradient of each component of the vector operator ${\widehat{\mathbf{\Theta }}}_{\mathbf{k}}$ applied to the Bloch state ${\mathbf{u}}_k$:(16) $$\begin{array}{cc}\hfill {\left(\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k\right)}^t·\mathbf{n}& =\left(\begin{array}{c}{n}_i\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\\ n_i\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\\ n_i\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\end{array}\right)\hfill \\ \hfill & =\left(\begin{array}{c}\mathbf{n}·{\mathbf{\nabla }}_{\mathbf{k}}{v}_x\\ \mathbf{n}·{\mathbf{\nabla }}_{\mathbf{k}}{v}_y\\ \mathbf{n}·{\mathbf{\nabla }}_{\mathbf{k}}{v}_z\end{array}\right)\text{,}\hfill \end{array}$$(16)

where Einstein’s notations are used and with $v_j\equiv {\mathbf{e}}_j·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k$.

Let us now determine the elements of the tensor $[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}]{\mathbf{u}}_k$ by calculating$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k& =\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}\left[(i\mathbf{k}+\mathrm{\nabla })\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k\right)\right]\hfill \\ \hfill & =\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}(i\mathbf{k}+\mathrm{\nabla })\times \left[\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k\right]\hfill \\ \hfill & +(i\mathbf{k}+\mathrm{\nabla })\times \frac{\mathit{\partial }}{\mathit{\partial }{k}_i}\left[\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k\right]\hfill \\ \hfill & =i{\mathbf{e}}_i\times \left[\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k\right]\hfill \\ \hfill & +(i\mathbf{k}+\mathrm{\nabla })\times \left[\frac{1}{\mathit{\epsilon }(\mathbf{r})}i{\mathbf{e}}_i\times {\mathbf{u}}_k\right]\text{.}\hfill \end{array}$$

Let us rewrite this result as$$\frac{\mathit{\partial }}{\mathit{\partial }{k}_j}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k=i{\mathbf{e}}_i\times \mathbf{b}+(i\mathbf{k}+\mathrm{\nabla })\times i{\mathbf{A}}_i\text{,}$$

with $\mathbf{b}=\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k$ and ${\mathbf{A}}_i=\frac{1}{\mathit{\epsilon }(\mathbf{r})}{\mathbf{e}}_i\times {\mathbf{u}}_k$.

By injecting the above expression in (15(15) $$\begin{array}{c}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k=\sum _{i,j}\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_i}{\mathbf{e}}_j·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}{\mathbf{u}}_k\right){\mathbf{e}}_i{\mathbf{e}}_j\text{,}\end{array}$$(15) ), we obtain$$\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k=\sum _{i,j}\left[i{\mathbf{e}}_i\times \mathbf{b}·{\mathbf{e}}_j+(i\mathbf{k}+\mathrm{\nabla })\times i{\mathbf{A}}_i·{\mathbf{e}}_j\right]{\mathbf{e}}_i{\mathbf{e}}_j.$$

If we apply the definition of the vectorial product between the unitary tensor and an arbitrary vector, i.e. $I\times \mathbf{b}={\sum }_{i,j}\left({\mathbf{e}}_i\times \mathbf{b}·{\mathbf{e}}_j\right){\mathbf{e}}_i{\mathbf{e}}_j$ (cf. Appendix Equation (A8(A8) $$\begin{array}{c}\hfill I\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{e}}_x\times \mathbf{b})}^t\\ {({\mathbf{e}}_y\times \mathbf{b})}^t\\ {({\mathbf{e}}_z\times \mathbf{b})}^t\end{array}\right)=\left(\begin{array}{ccc}0& -b_z& b_y\\ b_z& 0& -b_x\\ -b_y& b_x& 0\end{array}\right)\text{.}\end{array}$$(A8) )), and the definition of the vectorial product between an arbitrary vector and an arbitrary tensor, i.e. $\mathbf{c}\times A={\sum }_{i,j}\left(\mathbf{c}\times {\mathbf{A}}_i·{\mathbf{e}}_j\right){\mathbf{e}}_i{\mathbf{e}}_j$ (cf. Appendix Equation (A7(A7) $$\begin{array}{c}\hfill \mathbf{b}\times A=\left(\begin{array}{c}{(\mathbf{b}\times {\mathbf{A}}_x)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_y)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_z)}^t\end{array}\right)\text{.}\end{array}$$(A7) )), we can rewrite the above expression in a compact form:(17) $$\begin{array}{cc}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k& =iI\times \mathbf{b}+(i\mathbf{k}+\mathrm{\nabla })\times iA\hfill \\ \hfill & =iI\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}(i\mathbf{k}+\mathrm{\nabla })\times {\mathbf{u}}_k\right)\hfill \\ \hfill & +(i\mathbf{k}+\mathrm{\nabla })\times \left(i\frac{1}{\mathit{\epsilon }(\mathbf{r})}I\times {\mathbf{u}}_k\right)\text{.}\hfill \end{array}$$(17)

Finally, the expression of the tensor operator defined by (14(14) $$\begin{array}{c}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]=\left(\begin{array}{ccc}\frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_x}\\ \frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_y}\\ \frac{\mathit{\partial }{\mathbf{e}}_x·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}& \frac{\mathit{\partial }{\mathbf{e}}_y·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}& \frac{\mathit{\partial }{\mathbf{e}}_z·{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}}{\mathit{\partial }{k}_z}\end{array}\right).\end{array}$$(14) ) is(18) $$\begin{array}{cc}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \left((i\mathbf{k}+\mathrm{\nabla })\times \right)\hfill \\ \hfill & +(i\mathbf{k}+\mathrm{\nabla })\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \right)\text{.}\hfill \end{array}$$(18)

It is important to note that the scalar field $\frac{1}{\mathit{\epsilon }(\mathbf{r})}$ can be put in front of the first term whereas it must stay after the $\mathrm{\nabla }\times$ operator in the second term.

6. Proof of the formula

In order to prove Equation (13(13) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})}{{\mathcal{U}}_H+{\mathcal{U}}_E}\text{.}\end{array}$$(13) ), let us first define the inner product between a vector (here ${\mathbf{u}}_k$) and a rank-two tensor (here $[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}]{\mathbf{u}}_k$) as (see also Notes)(19) $$\begin{array}{c}\hfill \left({\mathbf{u}}_k,\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k\right)\equiv \int d^3\mathbf{r}\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }\text{,}\end{array}$$(19)

where the dot symbol denotes the matrix product.

Let us, then, separate the tensorial operator (18(18) $$\begin{array}{cc}\hfill \left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \left((i\mathbf{k}+\mathrm{\nabla })\times \right)\hfill \\ \hfill & +(i\mathbf{k}+\mathrm{\nabla })\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \right)\text{.}\hfill \end{array}$$(18) ) into two terms:(20) $$\begin{array}{c}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1=\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[iI\times \left(i\mathbf{k}\times \right)+i\mathbf{k}\times \left(iI\times \right)\right]\end{array}$$(20)

and(21) $$\begin{array}{c}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2=\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times )+\mathrm{\nabla }\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \right)\text{.}\end{array}$$(21)

The first term is a tensor operator involving $\mathbf{k}$ without spatial derivatives (for this reason, $\frac{1}{\mathit{\epsilon }(\mathbf{r})}$ can be put in front of the operator) whereas the second term is a pure differential tensor operator. We note that ${\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1^{\ast }={\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1$ whereas ${\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2^{\ast }=-{\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2$, provided $\mathit{\epsilon }$ is real (loss-less medium).

Let us calculate the inner product that is associated with the first operator, Equation (20(20) $$\begin{array}{c}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1=\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[iI\times \left(i\mathbf{k}\times \right)+i\mathbf{k}\times \left(iI\times \right)\right]\end{array}$$(20) ):(22) $$\begin{array}{cc}& {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[iI\times (i\mathbf{k}\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }+i\mathbf{k}\times (iI\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }\right]\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[i(i\mathbf{k}\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }+i(i\mathbf{k}\times {\mathbf{u}}_k^{\ast })\times {\mathbf{u}}_k\right]\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[i{\mathbf{u}}_k(i\mathbf{k}·{\mathbf{u}}_k^{\ast })+\mathbf{k}|{\mathbf{u}}_k{|}^2+i{\mathbf{u}}_k^{\ast }(i\mathbf{k}·{\mathbf{u}}_k)+\mathbf{k}{|{\mathbf{u}}_k|}^2\right]\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}\left[2\mathbf{k}|{\mathbf{u}}_k{|}^2+2\mathfrak{R}\left(i{\mathbf{u}}_k(\mathrm{\nabla }·{\mathbf{u}}_k^{\ast })\right)\right]\text{,}\hfill \end{array}$$(22)

where calculation rules (A10(A10) $$\begin{array}{c}\hfill I\times (\mathbf{b}\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{c})\times \mathbf{a}\text{.}\end{array}$$(A10) ) and (A11(A11) $$\begin{array}{c}\hfill \mathbf{b}\times (I\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{a})\times \mathbf{c}\text{.}\end{array}$$(A11) ) (see Appendix 1) were applied to, respectively, the first and second terms on the right-hand side and the transversality condition (6(6) $$\begin{array}{c}\hfill \mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})=-i\mathbf{k}·{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})\text{.}\end{array}$$(6) ) was used. The result is a real vector, as required. The contribution of the first part of the operator to the inner product is therefore:(23) $$\begin{array}{cc}& \left({\mathbf{u}}_k,{\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1{\mathbf{u}}_k\right)\hfill \\ \hfill & =\int d^3\mathbf{r}\frac{1}{\mathit{\epsilon }(\mathbf{r})}2\left[\mathbf{k}|{\mathbf{u}}_k{|}^2+\mathfrak{R}\left(i{\mathbf{u}}_k(\mathrm{\nabla }·{\mathbf{u}}_k^{\ast })\right)\right]\text{.}\hfill \end{array}$$(23)

Let us now calculate the inner product that is associated with the second operator, Equation (21(21) $$\begin{array}{c}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2=\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times )+\mathrm{\nabla }\times \left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times \right)\text{.}\end{array}$$(21) ):(24) $$\begin{array}{cc}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & +i\mathrm{\nabla }\times \left(I\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\right)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\hfill \\ \hfill & -i\left({\mathbf{u}}_k^{\ast }\times \mathrm{\nabla }\right)\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\text{,}\hfill \end{array}$$(24)

where calculation rules (A12(A12) $$\begin{array}{c}\hfill I\times (\mathrm{\nabla }\times \mathbf{c})·\mathbf{a}=(\mathrm{\nabla }\times \mathbf{c})\times \mathbf{a}\text{.}\end{array}$$(A12) ) and (A14(A14) $$\begin{array}{c}\hfill \mathrm{\nabla }\times (I\times \mathbf{c})·\mathbf{a}=-(\mathbf{a}\times \mathrm{\nabla })\times \mathbf{c}\text{.}\end{array}$$(A14) ) (see Appendix 1) were applied to, respectively, the first and second terms on the right-hand side.

Although nothing allows us to guess it, the second term on the right-hand side of the second equality in (24(24) $$\begin{array}{cc}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & +i\mathrm{\nabla }\times \left(I\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\right)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\hfill \\ \hfill & -i\left({\mathbf{u}}_k^{\ast }\times \mathrm{\nabla }\right)\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\text{,}\hfill \end{array}$$(24) ) is, after integration over the crystal unit cell, the complex conjugate of the first term, as we will show immediately.

Leaving on the side, but without forgetting it, the factor $-i$, let us develop the term $\left({\mathbf{u}}^{\ast }\times \mathrm{\nabla }\right)\times \frac{\mathbf{u}}{\mathit{\epsilon }}$ (subscript $\mathbf{k}$ and dependence on $\mathbf{r}$ are omitted for conciseness) by applying usual calculation rules of linear algebra (see Appendix 1).$$\begin{array}{cc}\hfill \left({\mathbf{u}}^{\ast }\times \mathrm{\nabla }\right)\times \frac{\mathbf{u}}{\mathit{\epsilon }}& =\left(\mathrm{\nabla }\frac{\mathbf{u}}{\mathit{\epsilon }}\right){\mathbf{u}}^{\ast }-\left(\mathrm{\nabla }·\frac{\mathbf{u}}{\mathit{\epsilon }}\right){\mathbf{u}}^{\ast }\hfill \\ \hfill & =\left[\mathrm{\nabla }\left({\mathbf{u}}^{\ast }·\frac{\mathbf{u}}{\mathit{\epsilon }}\right)-(\mathrm{\nabla }{\mathbf{u}}^{\ast })\frac{\mathbf{u}}{\mathit{\epsilon }}\right]\hfill \\ \hfill & -\left[\mathrm{\nabla }·\left({\mathbf{u}}^{\ast }\frac{\mathbf{u}}{\mathit{\epsilon }}\right)-\left(\frac{\mathbf{u}}{\mathit{\epsilon }}·\mathrm{\nabla }\right){\mathbf{u}}^{\ast }\right]\hfill \\ \hfill & =\mathrm{\nabla }\frac{{|\mathbf{u}|}^2}{\mathit{\epsilon }}-\mathrm{\nabla }·\left(\frac{{\mathbf{u}}^{\ast }\mathbf{u}}{\mathit{\epsilon }}\right)\hfill \\ \hfill & +\left(\frac{\mathbf{u}}{\mathit{\epsilon }}·\mathrm{\nabla }\right){\mathbf{u}}^{\ast }-(\mathrm{\nabla }{\mathbf{u}}^{\ast })\frac{\mathbf{u}}{\mathit{\epsilon }}\text{,}\hfill \end{array}$$

where (A15(A15) $$\begin{array}{c}\hfill (\mathbf{u}\times \mathrm{\nabla })\times \mathbf{v}=(\mathrm{\nabla }\mathbf{v})\mathbf{u}-(\mathrm{\nabla }·\mathbf{v})\mathbf{u}.\end{array}$$(A15) ) was used to write the first equality, (A18(A18) $$\begin{array}{c}\hfill (\mathrm{\nabla }\mathbf{v})\mathbf{u}=\mathrm{\nabla }(\mathbf{u}·\mathbf{v})-(\mathrm{\nabla }\mathbf{u})\mathbf{v}\text{.}\end{array}$$(A18) ) and (A17(A17) $$\begin{array}{c}\hfill (\mathrm{\nabla }·\mathbf{v})\mathbf{u}=\mathrm{\nabla }·(\mathbf{u}\mathbf{v})-(\mathbf{v}·\mathrm{\nabla })\mathbf{u}\text{.}\end{array}$$(A17) ) were used to develop, respectively, the first and second terms, and terms were rearranged to lead the last equality.

After integration over the unit cell volume $\mathrm{\Omega }$, the first two terms vanish individually. Indeed, by transforming integrals over the unit cell volume into integrals on the unit cell surface S, one obtains$$\begin{array}{cc}& {\int }_{\mathrm{\Omega }}\mathrm{\nabla }\frac{{|\mathbf{u}|}^2}{\mathit{\epsilon }}{d}^3\mathbf{r}-{\int }_{\mathrm{\Omega }}\mathrm{\nabla }·\left(\frac{{\mathbf{u}}^{\ast }\mathbf{u}}{\mathit{\epsilon }}\right)d^3\mathbf{r}\hfill \\ \hfill & ={\int }_S\frac{{|\mathbf{u}|}^2}{\mathit{\epsilon }}\mathbf{n}\mathrm{d}S-{\int }_S\left(\frac{{\mathbf{u}}^{\ast }\mathbf{u}}{\mathit{\epsilon }}\right)·\mathbf{n}\mathrm{d}S=0\hfill \end{array}$$

by virtue of discrete translational symmetry of ${\mathbf{u}}_{\mathbf{k}}(\mathbf{r}),{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})$ and $\mathit{\epsilon }(\mathbf{r})$ imposed by periodic boundary conditions on the unit cell surface.

The remaining terms can be regrouped using (A16(A16) $$\begin{array}{c}\hfill (\mathrm{\nabla }\times \mathbf{v})\times \mathbf{u}=(\mathbf{u}·\mathrm{\nabla })\mathbf{v}-(\mathrm{\nabla }\mathbf{v})\mathbf{u}\text{.}\end{array}$$(A16) ), which gives$$\left(\frac{\mathbf{u}}{\mathit{\epsilon }}·\mathrm{\nabla }\right){\mathbf{u}}^{\ast }-(\mathrm{\nabla }{\mathbf{u}}^{\ast })\frac{\mathbf{u}}{\mathit{\epsilon }}=(\mathrm{\nabla }\times {\mathbf{u}}^{\ast })\times \frac{\mathbf{u}}{\mathit{\epsilon }}=\frac{1}{\mathit{\epsilon }}(\mathrm{\nabla }\times {\mathbf{u}}^{\ast })\times \mathbf{u}$$

and, after multiplication by $-i$ (the factor left in standby above) and return to the Bloch state notation,$$-i\frac{1}{\mathit{\epsilon }}(\mathrm{\nabla }\times {\mathbf{u}}_k^{\ast })\times {\mathbf{u}}_k={\left(i\frac{1}{\mathit{\epsilon }}(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\right)}^{\ast }\text{,}$$

that is, the complex conjugate of the first term of the right-hand side of the second equality in (24(24) $$\begin{array}{cc}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & +i\mathrm{\nabla }\times \left(I\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\right)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\hfill \\ \hfill & -i\left({\mathbf{u}}_k^{\ast }\times \mathrm{\nabla }\right)\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\text{,}\hfill \end{array}$$(24) ).

Since both terms on the right-hand side of the second equality in (24(24) $$\begin{array}{cc}\hfill {\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k·{\mathbf{u}}_k^{\ast }& =\frac{1}{\mathit{\epsilon }(\mathbf{r})}iI\times (\mathrm{\nabla }\times {\mathbf{u}}_k)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & +i\mathrm{\nabla }\times \left(I\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\right)·{\mathbf{u}}_k^{\ast }\hfill \\ \hfill & =\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\hfill \\ \hfill & -i\left({\mathbf{u}}_k^{\ast }\times \mathrm{\nabla }\right)\times \frac{{\mathbf{u}}_k}{\mathit{\epsilon }(\mathbf{r})}\text{,}\hfill \end{array}$$(24) ) are complex conjugate of each others, the contribution of the second part of the operator to the inner product is(25) $$\begin{array}{cc}& \left({\mathbf{u}}_k,{\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k\right)\hfill \\ \hfill & =\int d^3\mathbf{r}2\mathfrak{R}\left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\right)\text{.}\hfill \end{array}$$(25)

The result is a real vector, as required.

If we now sum up the contributions of both operators, Equations (23(23) $$\begin{array}{cc}& \left({\mathbf{u}}_k,{\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_1{\mathbf{u}}_k\right)\hfill \\ \hfill & =\int d^3\mathbf{r}\frac{1}{\mathit{\epsilon }(\mathbf{r})}2\left[\mathbf{k}|{\mathbf{u}}_k{|}^2+\mathfrak{R}\left(i{\mathbf{u}}_k(\mathrm{\nabla }·{\mathbf{u}}_k^{\ast })\right)\right]\text{.}\hfill \end{array}$$(23) ) and (25(25) $$\begin{array}{cc}& \left({\mathbf{u}}_k,{\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]}_2{\mathbf{u}}_k\right)\hfill \\ \hfill & =\int d^3\mathbf{r}2\mathfrak{R}\left(\frac{1}{\mathit{\epsilon }(\mathbf{r})}i(\mathrm{\nabla }\times {\mathbf{u}}_k)\times {\mathbf{u}}_k^{\ast }\right)\text{.}\hfill \end{array}$$(25) ), and multiply by $\frac{c^2}{2\mathit{\omega }}$, we obtain an expression that is equal to the integral of the Poynting vector, Equation (11(11) $$\begin{array}{cc}\hfill \mathbf{S}(\mathbf{r})& =\frac{1}{2\mathit{\omega }{\mathit{\epsilon }}_0\mathit{\epsilon }(\mathbf{r})}\left(\mathbf{k}|{\mathbf{u}}_{\mathbf{k}}{(\mathbf{r})|}^2\right.\hfill \\ \hfill & +\mathfrak{R}\left[i{\mathbf{u}}_{\mathbf{k}}(\mathbf{r})(\mathrm{\nabla }·{\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r}))\right.\hfill \\ \hfill & \left.\left.+i(\mathrm{\nabla }\times {\mathbf{u}}_{\mathbf{k}}(\mathbf{r}))\times {\mathbf{u}}_{\mathbf{k}}^{\ast }(\mathbf{r})\right]\right)\text{.}\hfill \end{array}$$(11) ), multiplied by $\frac{2}{{\mathit{\mu }}_0}$.

This completes the proof of$$\frac{c^2}{2\mathit{\omega }}\left({\mathbf{u}}_k,\left[{\mathbf{\nabla }}_{\mathbf{k}}{\widehat{\mathbf{\Theta }}}_{\mathbf{k}}\right]{\mathbf{u}}_k\right)=\frac{2}{{\mathit{\mu }}_0}\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})$$

from which Equation (13(13) $$\begin{array}{c}\hfill {\mathbf{v}}_g=\frac{\int d^3\mathbf{r}\mathbf{S}(\mathbf{r})}{{\mathcal{U}}_H+{\mathcal{U}}_E}\text{.}\end{array}$$(13) ), our central result, follows.

7. Conclusion

We provided the proof of a less commonly used formula of the group velocity of Bloch modes in periodic media. The proof relies on the explicit definition and the obtained expression of the tensor operator resulting from differentiation with respect to wave-vector of the vector operator acting on Bloch modes. In addition to have immediate physical interpretation in terms of energy flow, this formula may facilitate group velocity calculation in photonics crystal materials and, in particular, those designed to exploit slow light.

8. Notes

For complex vector fields $\mathbf{F}(\mathbf{r})$ and $\mathbf{G}(\mathbf{r})$, the inner product (scalar quantity) is defined by $(\mathbf{F},\mathbf{G})=\int d^3\mathbf{r}{\mathbf{F}}^{\ast }(\mathbf{r})·\mathbf{G}(\mathbf{r})$, the volume integral being performed on the crystal unit cell due to periodic boundary conditions. This definition of the inner product between two vector fields is easily generalized to the case involving a vector field $\mathbf{F}(\mathbf{r})$ and a second-rank tensor field $\mathbf{T}(\mathbf{r})$, the result being this time a vector quantity. In a matrix representation of vectors and tensors in three-dimensional space, the inner product takes the following form: $(\mathbf{F},\mathbf{T})=\int d^3\mathbf{r}\mathbf{T}(\mathbf{r})·{\mathbf{F}}^{\ast }(\mathbf{r})$, where the dot symbol denotes the matrix product. In the particular case where $\mathbf{T}(\mathbf{r})=\mathbf{v}\mathbf{F}(\mathbf{r})$ with $\mathbf{v}$ a constant vector, we have $(\mathbf{F},\mathbf{v}F)=\int d^3{\mathbf{r}\mathbf{v}|\mathbf{F}|}^2=\mathbf{v}\int d^3\mathbf{r}{|\mathbf{F}|}^2=\mathbf{v}(\mathbf{F},\mathbf{F})$.

Acknowledgements

Olivier Deparis thanks Alexandre Mauroy (Department of Mathematics, University of Namur) for enlightening discussion related to algebraic aspects of the present study.

Notes

No potential conflict of interest was reported by the authors.

Appendix 1

Let us note $\left\{{\mathbf{e}}_x,{\mathbf{e}}_y,{\mathbf{e}}_z\right\}$ the cartesian vector basis in which all the following calculation rules are established. In this basis, let us express an arbitrary vector by $\mathbf{b}={\sum }_k{b}_k{\mathbf{e}}_k$. In a matrix form, we choose to represent $\mathbf{b}$ as a column vector, i.e. a $3\times 1$ matrix. Let $A={\sum }_{i,j}{A}_{ij}{\mathbf{e}}_i{\mathbf{e}}_j$ be an arbitrary rank-two tensor. The dyad ${\mathbf{e}}_i{\mathbf{e}}_j$ can be represented by a $3\times 3$ matrix with all matrix elements equal to zero, except one matrix element that is equal to unity. By convention, this unitary element is placed at the intersection of row i and column j of the matrix.

We begin by establishing an operational definition of $A\times \mathbf{b}$: that is, we aim to generalize the vector (cross) product ($\times$) to an operation between a rank-two tensor and a rank-one tensor (i.e. vector). Using Einstein’s notation for the indices $\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\in x,y,z$, the resulting rank-two tensor can be written as:(A1) $$\begin{array}{c}\hfill A\times \mathbf{b}=A_{\mathit{\alpha }\mathit{\beta }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\times b_{\mathit{\gamma }}{\mathbf{e}}_{\mathit{\gamma }}=A_{\mathit{\alpha }\mathit{\beta }}{b}_{\mathit{\gamma }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\times {\mathbf{e}}_{\mathit{\gamma }}\text{.}\end{array}$$(A1)

The vector product between two basis vectors can be expressed conveniently using Levi–Civita notation:(A2) $$\begin{array}{c}\hfill {\mathbf{e}}_{\mathit{\beta }}\times {\mathbf{e}}_{\mathit{\gamma }}={\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\mu }}{\mathbf{e}}_{\mathit{\mu }}\text{,}\end{array}$$(A2)

where ${\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\mu }}$ is equal to $0,+1,-1$ if $\mathit{\beta }=\mathit{\gamma }$, $\mathit{\beta }<\mathit{\gamma }$, $\mathit{\beta }>\mathit{\gamma }$, respectively, in the circular permutation of the indices. Using this notation, (A2(A2) $$\begin{array}{c}\hfill {\mathbf{e}}_{\mathit{\beta }}\times {\mathbf{e}}_{\mathit{\gamma }}={\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\mu }}{\mathbf{e}}_{\mathit{\mu }}\text{,}\end{array}$$(A2) ) becomes(A3) $$\begin{array}{c}\hfill A\times \mathbf{b}=A_{\mathit{\alpha }\mathit{\beta }}{b}_{\mathit{\gamma }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\times {\mathbf{e}}_{\mathit{\gamma }}=A_{\mathit{\alpha }\mathit{\gamma }}{b}_{\mathit{\delta }}{\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\delta }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\text{.}\end{array}$$(A3)

The second equality results from circular permutation $(\mathit{\beta },\mathit{\gamma },\mathit{\mu })\to (\mathit{\mu },\mathit{\beta },\mathit{\gamma })$ and substitution $\mathit{\beta }\to \mathit{\gamma }$, $\mathit{\mu }\to \mathit{\beta }$, $\mathit{\gamma }\to \mathit{\delta }$. The explicit matrix representation of (A3(A3) $$\begin{array}{c}\hfill A\times \mathbf{b}=A_{\mathit{\alpha }\mathit{\beta }}{b}_{\mathit{\gamma }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\times {\mathbf{e}}_{\mathit{\gamma }}=A_{\mathit{\alpha }\mathit{\gamma }}{b}_{\mathit{\delta }}{\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\delta }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}\text{.}\end{array}$$(A3) ) is therefore:(A4) $$\begin{array}{c}\hfill A\times \mathbf{b}=\left(\begin{array}{ccc}{A}_{xy}{b}_z-A_{xz}{b}_y& A_{xz}{b}_x-A_{xx}{b}_z& A_{xx}{b}_y-A_{xy}{b}_x\\ A_{yy}{b}_z-A_{yz}{b}_y& A_{yz}{b}_x-A_{yx}{b}_z& A_{yx}{b}_y-A_{yy}{b}_x\\ A_{zy}{b}_z-A_{zz}{b}_y& A_{zz}{b}_x-A_{zx}{b}_z& A_{zx}{b}_y-A_{zy}{b}_x\text{.}\end{array}\right).\end{array}$$(A4)

Relying on this matrix representation, we introduce an operational definition of $A\times \mathbf{b}$ using three vectors extracted from the matrix: ${\mathbf{A}}_j^t\equiv (A_{jx},A_{jy},A_{jz})$ where $j=x,y,z$ (superscript t denotes matrix transpose operation). In an explicit manner, we represent the rank-two tensor as(A5) $$\begin{array}{c}\hfill A=\left(\begin{array}{ccc}{A}_{xx}& A_{xy}& A_{xz}\\ A_{yx}& A_{yy}& A_{yz}\\ A_{zx}& A_{zy}& A_{zz}\end{array}\right)=\left(\begin{array}{c}{\mathbf{A}}_x^t\\ {\mathbf{A}}_y^t\\ {\mathbf{A}}_z^t\end{array}\right)\text{.}\end{array}$$(A5)

With these notations in hand, using usual vector and matrix calculus rules, it is straightforward to check that (A4(A4) $$\begin{array}{c}\hfill A\times \mathbf{b}=\left(\begin{array}{ccc}{A}_{xy}{b}_z-A_{xz}{b}_y& A_{xz}{b}_x-A_{xx}{b}_z& A_{xx}{b}_y-A_{xy}{b}_x\\ A_{yy}{b}_z-A_{yz}{b}_y& A_{yz}{b}_x-A_{yx}{b}_z& A_{yx}{b}_y-A_{yy}{b}_x\\ A_{zy}{b}_z-A_{zz}{b}_y& A_{zz}{b}_x-A_{zx}{b}_z& A_{zx}{b}_y-A_{zy}{b}_x\text{.}\end{array}\right).\end{array}$$(A4) ) can be rewritten as:(A6) $$\begin{array}{c}\hfill A\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{A}}_x\times \mathbf{b})}^t\\ {({\mathbf{A}}_y\times \mathbf{b})}^t\\ {({\mathbf{A}}_z\times \mathbf{b})}^t\end{array}\right)\text{.}\end{array}$$(A6)

This is just the operational definition we searched for.

An operational definition of $\mathbf{b}\times A$ can be established following the same arguments:(A7) $$\begin{array}{c}\hfill \mathbf{b}\times A=\left(\begin{array}{c}{(\mathbf{b}\times {\mathbf{A}}_x)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_y)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_z)}^t\end{array}\right)\text{.}\end{array}$$(A7)

Comparing (A6(A6) $$\begin{array}{c}\hfill A\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{A}}_x\times \mathbf{b})}^t\\ {({\mathbf{A}}_y\times \mathbf{b})}^t\\ {({\mathbf{A}}_z\times \mathbf{b})}^t\end{array}\right)\text{.}\end{array}$$(A6) ) and (A7(A7) $$\begin{array}{c}\hfill \mathbf{b}\times A=\left(\begin{array}{c}{(\mathbf{b}\times {\mathbf{A}}_x)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_y)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_z)}^t\end{array}\right)\text{.}\end{array}$$(A7) ), we find immediately that $\mathbf{b}\times A=-A\times \mathbf{b}$. This result could be inferred from our generalization of the vector product, as well. We have to be careful, however, that $({\mathbf{e}}_{\mathit{\alpha }}\times {\mathbf{e}}_{\mathit{\beta }}){\mathbf{e}}_{\mathit{\gamma }}={\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\mu }}^{\prime }{\mathbf{e}}_{\mathit{\mu }}{\mathbf{e}}_{\mathit{\alpha }}$ where ${\mathit{\epsilon }}_{\mathit{\beta }\mathit{\gamma }\mathit{\mu }}^{\prime }$ is here equal to $0,-1,+1$ if $\mathit{\beta }=\mathit{\gamma }$, $\mathit{\beta }<\mathit{\gamma }$, $\mathit{\beta }>\mathit{\gamma }$, respectively. This difference is responsible for the change of sign in $\mathbf{b}\times A=-A\times \mathbf{b}$.

In the particular case where the rank-two tensor is the identity tensor $\mathbf{I}={\mathit{\delta }}_{\mathit{\alpha }\mathit{\beta }}{\mathbf{e}}_{\mathit{\alpha }}{\mathbf{e}}_{\mathit{\beta }}$ (${\mathit{\delta }}_{\mathit{\alpha }\mathit{\beta }}$: Kronecker symbol), applying the operational definition (A6(A6) $$\begin{array}{c}\hfill A\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{A}}_x\times \mathbf{b})}^t\\ {({\mathbf{A}}_y\times \mathbf{b})}^t\\ {({\mathbf{A}}_z\times \mathbf{b})}^t\end{array}\right)\text{.}\end{array}$$(A6) ) leads to(A8) $$\begin{array}{c}\hfill I\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{e}}_x\times \mathbf{b})}^t\\ {({\mathbf{e}}_y\times \mathbf{b})}^t\\ {({\mathbf{e}}_z\times \mathbf{b})}^t\end{array}\right)=\left(\begin{array}{ccc}0& -b_z& b_y\\ b_z& 0& -b_x\\ -b_y& b_x& 0\end{array}\right)\text{.}\end{array}$$(A8)

We note that $I\times \mathbf{b}$ is an antisymmetric tensor. Of course, we have also $\mathbf{b}\times I=-I\times \mathbf{b}$.

Using matrix calculus rules and the above operational definitions, it is then easy to derive the following relations in which the dot symbol denotes the matrix product.

From (A8(A8) $$\begin{array}{c}\hfill I\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{e}}_x\times \mathbf{b})}^t\\ {({\mathbf{e}}_y\times \mathbf{b})}^t\\ {({\mathbf{e}}_z\times \mathbf{b})}^t\end{array}\right)=\left(\begin{array}{ccc}0& -b_z& b_y\\ b_z& 0& -b_x\\ -b_y& b_x& 0\end{array}\right)\text{.}\end{array}$$(A8) ), we deduce immediately(A9) $$\begin{array}{c}\hfill I\times \mathbf{b}·\mathbf{a}=\mathbf{b}\times \mathbf{a}\text{.}\end{array}$$(A9)

Substituting $\mathbf{b}\times \mathbf{c}$ to $\mathbf{b}$ in (A9(A9) $$\begin{array}{c}\hfill I\times \mathbf{b}·\mathbf{a}=\mathbf{b}\times \mathbf{a}\text{.}\end{array}$$(A9) ), we obtain another useful relation:(A10) $$\begin{array}{c}\hfill I\times (\mathbf{b}\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{c})\times \mathbf{a}\text{.}\end{array}$$(A10)

Using (A8(A8) $$\begin{array}{c}\hfill I\times \mathbf{b}=\left(\begin{array}{c}{({\mathbf{e}}_x\times \mathbf{b})}^t\\ {({\mathbf{e}}_y\times \mathbf{b})}^t\\ {({\mathbf{e}}_z\times \mathbf{b})}^t\end{array}\right)=\left(\begin{array}{ccc}0& -b_z& b_y\\ b_z& 0& -b_x\\ -b_y& b_x& 0\end{array}\right)\text{.}\end{array}$$(A8) ) to find the matrix representation of $I\times \mathbf{c}$, applying (A7(A7) $$\begin{array}{c}\hfill \mathbf{b}\times A=\left(\begin{array}{c}{(\mathbf{b}\times {\mathbf{A}}_x)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_y)}^t\\ {(\mathbf{b}\times {\mathbf{A}}_z)}^t\end{array}\right)\text{.}\end{array}$$(A7) ) to this result with $A\equiv I\times \mathbf{c}$ and multiplying the obtained matrix by a vector $\mathbf{a}$, we find(A11) $$\begin{array}{c}\hfill \mathbf{b}\times (I\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{a})\times \mathbf{c}\text{.}\end{array}$$(A11)

Similar relations can be derived when vector $\mathbf{b}$ is replaced by nabla differential operator $\mathrm{\nabla }\equiv \frac{\mathit{\partial }}{\mathit{\partial }x}{\mathbf{e}}_x+\frac{\mathit{\partial }}{\mathit{\partial }y}{\mathbf{e}}_y+\frac{\mathit{\partial }}{\mathit{\partial }z}{\mathbf{e}}_z$, provided that care is taken in order to respect the order involved in differential operations.

If $\mathrm{\nabla }$ is substituted for $\mathbf{b}$ in (A10(A10) $$\begin{array}{c}\hfill I\times (\mathbf{b}\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{c})\times \mathbf{a}\text{.}\end{array}$$(A10) ), the differential operator acts on $\mathbf{c}$ in a consistent manner on both the left-hand and right-hand sides of the equality. Therefore, we can safely proceed to the substitution and deduce the following relation:(A12) $$\begin{array}{c}\hfill I\times (\mathrm{\nabla }\times \mathbf{c})·\mathbf{a}=(\mathrm{\nabla }\times \mathbf{c})\times \mathbf{a}\text{.}\end{array}$$(A12)

If we intend to perform the same substitution from (A11(A11) $$\begin{array}{c}\hfill \mathbf{b}\times (I\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{a})\times \mathbf{c}\text{.}\end{array}$$(A11) ), we must first rewrite this relation by rearranging its right-hand side in such a way that, if $\mathrm{\nabla }$ is substituted for $\mathbf{b}$, it acts on $\mathbf{c}$ (and not on $\mathbf{a}$) consistently on both sides:(A13) $$\begin{array}{c}\hfill \mathbf{b}\times (I\times \mathbf{c})·\mathbf{a}=(\mathbf{b}\times \mathbf{a})\times \mathbf{c}=-(\mathbf{a}\times \mathbf{b})\times \mathbf{c}\text{.}\end{array}$$(A13)

This precaution being taken, we can safely proceed to the substitution and deduce the following relation:(A14) $$\begin{array}{c}\hfill \mathrm{\nabla }\times (I\times \mathbf{c})·\mathbf{a}=-(\mathbf{a}\times \mathrm{\nabla })\times \mathbf{c}\text{.}\end{array}$$(A14)

Other useful relations are listed hereafter without further comments since they can be either easily verified or found in textbookson linear algebra (for clarity, the dot symbol denotes exclusively the scalar product and the matrix product is not indicated by any symbol).(A15) $$\begin{array}{c}\hfill (\mathbf{u}\times \mathrm{\nabla })\times \mathbf{v}=(\mathrm{\nabla }\mathbf{v})\mathbf{u}-(\mathrm{\nabla }·\mathbf{v})\mathbf{u}.\end{array}$$(A15) (A16) $$\begin{array}{c}\hfill (\mathrm{\nabla }\times \mathbf{v})\times \mathbf{u}=(\mathbf{u}·\mathrm{\nabla })\mathbf{v}-(\mathrm{\nabla }\mathbf{v})\mathbf{u}\text{.}\end{array}$$(A16) (A17) $$\begin{array}{c}\hfill (\mathrm{\nabla }·\mathbf{v})\mathbf{u}=\mathrm{\nabla }·(\mathbf{u}\mathbf{v})-(\mathbf{v}·\mathrm{\nabla })\mathbf{u}\text{.}\end{array}$$(A17) (A18) $$\begin{array}{c}\hfill (\mathrm{\nabla }\mathbf{v})\mathbf{u}=\mathrm{\nabla }(\mathbf{u}·\mathbf{v})-(\mathrm{\nabla }\mathbf{u})\mathbf{v}\text{.}\end{array}$$(A18)