text:
# Canonical quantization of systems with time-dependent constraints11footnote 1e-mail: 𝑠𝑎𝑚𝑖₋𝑚𝑢𝑠𝑙𝑖ℎ@hotmail.com

## 1 Introduction

There are many singular theories with constraints which depend on time manifestly. The relativistic particle theories, string theories and theories of gravity are some examples. The presence of first class constraints in these theories forces one to impose time dependent gauge fixing condition for each first class constraint such that the whole set of constraints become converted in to second-class constraints. The rules of canonical quantization of systems with second-class constraints , which depend on time should be modified . However, there exist cases where one can make canonical transformation to make these constraints time-independent and then use the ordinary rules of quantization.

Recently, the canonical method \[4-7\] has been developed to investigate constrained systems. The equations of motion are obtained as total differential equations in many variables which require the investigation of integrability conditions. If the system is integrable, one can solve the equations of motion without using any gauge fixing conditions.

In this paper we consider the quantization of the relativistic spinless particle in the external field of a plane wave using the Hamilton-Jacobi method \[8-10\]. in fact, this work is a continuation of a previous work , where we have obtained the path integral for a relativistic charged particle in an external electromagnetic field and it is shown that the problems which arise from identifying the measure of the path integral when applying methods  and , are solved naturally by the canonical path integral formulation \[8-10\].

Now we would like to give a brief discussion of the canonical method. This method gives the set of Hamilton - Jacobi partial differential equations \[HJPDE\] as

$`H_\alpha ^{^{}}(t_\beta ,q_a,{\displaystyle \frac{S}{q_a}},{\displaystyle \frac{S}{t_a}})=0,`$
$`\alpha ,\beta =0,nr+1,\mathrm{},n,a=1,\mathrm{},nr,`$ (1)
where

$$H_\alpha ^{^{}}=H_\alpha (t_\beta ,q_a,p_a)+p_\alpha ,$$
(2)
and $`H_0`$ is defined as

$`H_0=p_aw_a+p_\mu \dot{q_\mu }|_{p_\nu =H_\nu }L(t,q_i,\dot{q_\nu },\dot{q_a}=w_a),`$
$`\mu ,\nu =nr+1,\mathrm{},n.`$ (3)
The equations of motion are obtained as total differential equations in many variables as follows:

$`dq_a={\displaystyle \frac{H_\alpha ^{^{}}}{p_a}}dt_\alpha ,dp_a={\displaystyle \frac{H_\alpha ^{^{}}}{q_a}}dt_\alpha ,dp_\beta ={\displaystyle \frac{H_\alpha ^{^{}}}{t_\beta }}dt_\alpha .`$ (4)
$`dz=(H_\alpha +p_a{\displaystyle \frac{H_\alpha ^{^{}}}{p_a}})dt_\alpha ;`$ (5)
$`\alpha ,\beta =0,nr+1,\mathrm{},n,a=1,\mathrm{},nr`$
where $`z=S(t_\alpha ;q_a)`$. The set of equations (4,5) is integrable  if

$$dH_0^{^{}}=0,dH_\mu ^{^{}}=0,\mu =np+1,\mathrm{},n.$$
(6)
If condition (6) are not satisfied identically, one considers them as new constraints and again testes the consistency conditions. Hence, the canonical formulation leads to obtain the set of canonical phase space coordinates $`q_a`$ and $`p_a`$ as functions of $`t_\alpha `$, besides the canonical action integral is obtained in terms of the canonical coordinates.The Hamiltonians $`H_\alpha ^{^{}}`$ are considered as the infinitesimal generators of canonical transformations given by parameters $`t_\alpha `$ respectively.

## 2 Quantization of constrained systems

For the quantization of constrained systems we can use the Dirac’s method of quantization , or the path integral quantization method \[8-10\].

Now will shall give a brief information about these two methods.

### 2.1 Operator quantization

For the Dirac’s quantization method we have

$$H_\alpha ^{^{}}\mathrm{\Psi }=0,\alpha =0,nr+1,\mathrm{},n,$$
(7)
where $`\mathrm{\Psi }`$ is the wave function. The consistency conditions are

$$[H_\mu ^{},H_\nu ^{}]\mathrm{\Psi }=0,\mu ,\nu =1,\mathrm{},r,$$
(8)
where$`[,]`$ is the commutator. The constraints $`H_\alpha ^{}`$ are called first- class constraints if they satisfy

$$[H_\mu ^{},H_\nu ^{}]=C_{\mu \nu }^\gamma H_\gamma ^{}.$$
(9)
In the case when the Hamiltonians $`H_\mu ^{}`$ satisfy

$$[H_\mu ^{},H_\nu ^{}]=C_{\mu \nu },$$
(10)
with $`C_{\mu \nu }`$ do not depend on $`q_i`$ and $`p_i`$, then from (8) there arise naturally Dirac’ brackets and the canonical quantization will be performed taking Dirac’s brackets into commutators.

### 2.2 Path integral quantization method

The path integral quantization is an alternative method to perform the quantization of constrained systems.

Now we shall give a brief review of the canonical path integral formulation of constrained systems \[8-10\].

If the set of equations (4) is integrable then one can solve them to obtain the canonical phase-space coordinates as

$$q_aq_a(t,t_\mu ),p_ap_a(t,t_\mu ),\mu =1,\mathrm{},r,$$
(11)
In this case, the path integral representation may be written as \[8-10\]

$`Out|S|In={\displaystyle \underset{a=1}{\overset{nr}{}}dq^adp^a\mathrm{exp}[i\{_{t_\alpha }^{t_{}^{}{}_{\alpha }{}^{}}(H_\alpha +p_a\frac{H_\alpha ^{^{}}}{p_a})𝑑t_\alpha \}]},`$
$`a=1,\mathrm{},nr,\alpha =0,nr+1,\mathrm{},n.`$ (12)
One should notice that the integral (12) is an integration over the canonical phase - space coordinates $`(q_a,p_a)`$.

## 3 Quantization of the relativistic particle in a plane wave

Motion of spinless particle in the external field of a plane wave is defined by a singular Lagrangians. The usual approach that authors investigate is the Dirac’s approach .

Let us consider the reparametrization invariant action of a spinless particle in the external electromagnatic field of a plane wave directed a long the axis $`x^3`$ as

$`S=`$ $`{\displaystyle (m\sqrt{\dot{x}^2}+e\dot{x}A)𝑑\tau },`$ (13)
$`=`$ $`{\displaystyle [m\sqrt{2\dot{x}_+\dot{x}_{}(\dot{x}_{})^2}+e\dot{x}^aA_a(x_{})]𝑑\tau },`$
$`A_\mu =`$ $`(0,A^a,0),x_\pm =(x^0\pm x^3)/\sqrt{2},x_{}=(x^a),a=1,2.`$
The momenta conjugated to $`x_\pm `$ and $`x^a`$ are

$`\pi _+={\displaystyle \frac{L}{\dot{x}_+}}={\displaystyle \frac{m\dot{x}_{}}{\sqrt{2\dot{x}_+\dot{x}_{}(\dot{x}_{})^2}}},`$ (14)
$`\pi _{}={\displaystyle \frac{L}{\dot{x}_{}}}={\displaystyle \frac{m\dot{x}_+}{\sqrt{2\dot{x}_+\dot{x}_{}(\dot{x}_{})^2}}},`$ (15)
$`\pi _a={\displaystyle \frac{L}{\dot{x}^a}}={\displaystyle \frac{m\dot{x}^a}{\sqrt{2\dot{x}_+\dot{x}_{}(\dot{x}_{})^2}}}eA_a(x_{}),`$ (16)
the primary constraint is in the form

$$\varphi ^{(1)}=\pi _{}\frac{[\pi _a+eA_a(x_{})]^2+m^2}{2\pi _+}=0,\pi _\pm 0.$$
(17)
The canonical Hamiltonian (3) vanishes identically. So the primary Hamiltonian is the constraint itself multiplied by a function $`\lambda (\tau )`$, i.e

$$H^{(1)}=\lambda \varphi ^{(1)}.$$
(18)
The equations of motion read as

$`\dot{x}_+=\{x_+,H^{(1)}\}=\lambda {\displaystyle \frac{[\pi _a+eA_a(x_{})]^2+m^2}{\pi _+}},`$ (19)
$`\dot{x}_{}=\{x_{},H^{(1)}\}=\lambda ,`$ (20)
$`\dot{x}_a=\{x_a,H^{(1)}\}=\lambda {\displaystyle \frac{(\pi _a+eA_a(x_{}))}{\pi _+}},`$ (21)
$`\dot{\pi }_+=\{\pi _+,H^{(1)}\}=0,`$ (22)
$`\dot{\pi }_a=\{\pi _a,H^{(1)}\}=0,`$ (23)
$`\dot{\pi }_{}=\{\pi _{},H^{(1)}\}=\lambda e{\displaystyle \frac{(\pi _a+eA_a(x_{}))}{\pi _+}}{\displaystyle \frac{A_a}{x_{}}},`$ (24)
Since we have only one first class constraint, the theory is degenerate and it is impossible to find $`\lambda `$ in the Dirac’s procedure. To determine $`\lambda `$, one should impose time-dependent gauge fixing of the form

$$\mathrm{\Phi }^G=x_{}\tau =0.$$
(25)
From the condition of conservation of this gauge in the time $`\tau `$

$$\frac{d\mathrm{\Phi }^G}{d\tau }=\frac{\mathrm{\Phi }^G}{\tau }+\{\mathrm{\Phi }^G,H^{(1)}\}=(\lambda 1)=0,$$
(26)
we obtain $`\lambda =1`$. No other constraints arise.

Now we would like to study the quantization of the same problem using the canonical method \[4-7\]. Equations (14) and (16) lead us to obtain the expressible velocities $`\dot{x}_+`$ and $`\dot{x}_a`$ in terms of the primary unexpressible velocity $`\dot{x}_{}`$ as

$$\dot{x}_+=\frac{\pi _{}}{\pi _+}\dot{x}_{}=w_+,\dot{x}_a=\frac{(\pi _a+eA_a)}{\pi _+}\dot{x}_{}=w_a.$$
(27)
Substituting (27) in (15), one gets

$$\pi _{}=\frac{[\pi _a+eA_a(x_{})]^2+m^2}{2\pi _+}=H_x_{},\pi _\pm 0.$$
(28)
The canonical Hamiltonian $`H_\tau `$ can be written as

$$H_\tau =L+\pi _aw_a+\pi _+w_+\dot{x}_{}H_x_{}.$$
(29)
Explicit calculations show that $`H_\tau `$ vanishes identically.

Making use of equations (2) and (28,29), one obtains the set of Hamilton Jacob partial differential equations as

$`H_{}^{}{}_{\tau }{}^{}=\pi _\tau =0,\pi _\tau ={\displaystyle \frac{S}{\tau }},`$ (30)
$`H_{}^{}{}_{x_{}}{}^{}=\pi _{}+H_x_{}=0,\pi _{}={\displaystyle \frac{S}{x_{}}}.`$ (31)
The equations of motion are obtained as total differential equations in many variables as follows:

$`dx_+`$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{\pi _+}}d\tau +{\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{\pi _+}}dx_{},`$ (32)
$`={\displaystyle \frac{[\pi _a+eA_a(x_{})]^2+m^2}{\pi _+}}dx_{},`$
$`dx_a`$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{\pi _a}}d\tau +{\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{\pi _a}}dx_{},`$ (33)
$`={\displaystyle \frac{(\pi _a+eA_a(x_{}))}{\pi _+}}dx_{},`$
$`d\pi _+`$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{x_+}}d\tau {\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{x_+}}dx_{}=0,`$ (34)
$`d\pi _a`$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{x_a}}d\tau {\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{x_a}}dx_{}=0,`$ (35)
$`d\pi _{}`$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{x_{}}}d\tau {\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{x_{}}}dx_{}=0,`$ (37)
$`=e{\displaystyle \frac{(\pi _a+eA_a(x_{}))}{\pi _+}}{\displaystyle \frac{A_a}{x_{}}}dx_{},`$
$`d\pi _\tau `$ $`={\displaystyle \frac{H_{}^{}{}_{\tau }{}^{}}{\tau }}d\tau {\displaystyle \frac{H_{}^{}{}_{x_{}}{}^{}}{\tau }}dx_{}=0.`$ (38)
In order to have a consistent theory, one should consider the variations of constraints in general. Since the total variations of the constraints $`H_{}^{}{}_{\tau }{}^{}`$ and $`H_{}^{}{}_{x_{}}{}^{}`$ are identically zero, no further constraints arise and the equations of motion are integrable. Hence, the canonical phase space coordinates $`(x_+,\pi _+)`$ and $`(x_a,\pi _a)`$ are obtained in terms of parameters $`\tau `$ and $`x_{}`$.

To obtain the operator quantization of this system one can follow the procedure discussed in section (2.1). In this case one takes the constraint equation as an operator whose action on the allowed Hilbert space vectors is constrained to $`zero`$, i., e., $`H_{}^{}{}_{x_{}}{}^{}\mathrm{\Psi }=0`$, we obtain

$$i\frac{\mathrm{\Psi }}{x_{}}=\frac{[i\frac{}{x^a}eA^a(x_{})]^2+m^2}{2\pi _+}\mathrm{\Psi }$$
(39)
One should notice that (39) is the Klein-Gordon equation. In fact, the Klein-Gordon equation

$$(𝒫^2m^2)\mathrm{\Psi }(x)=0,𝒫_\mu =i_\mu eA_\mu (x),$$
(40)
in the external field of a plane wave (13), being written in the light cone variables $`(x_\pm ,x_{}=(x^a))`$, reads as

$$(2\frac{^2}{x_{}x_+}+[i\frac{}{x^a}eA^a(x_{})]^2+m^2)\mathrm{\Psi }(x)=0.$$
(41)
Now to obtain the path integral quantization of this system, we can use equation (5) to obtain the canonical action as

$$S=(\frac{[\pi _a+eA_a(x_{})]^2+m^2}{2\pi _+}+\pi _+\dot{x}_++\pi _a\dot{x}_a)𝑑x_{}.$$
(42)
Making use of (42) and (12) the path integral for the system (13) is obtained as

$`x_+,x_a,\tau ,x_{};x_{}^{}{}_{+}{}^{},x_{}^{}{}_{a}{}^{},\tau ^{^{}},x_{}^{}{}_{}{}^{}={\displaystyle _{x_+,x^a}^{x_{}^{}{}_{+}{}^{},x_{}^{}{}_{a}{}^{}}}{\displaystyle dx_+dx_ad\pi _+d\pi _a}`$
$`\mathrm{exp}[i\{{\displaystyle _x_{}^x_{}^{}}({\displaystyle \frac{[\pi _a+eA_a(x_{})]^2+m^2}{2\pi _+}}+\pi _+\dot{x}_++\pi _a\dot{x}_a)𝑑x_{}\}].`$ (43)
In fact, the path integral representation (43) is an integration over the canonical phase-space coordinates This path integral representation is an integration over the canonical phase space coordinates $`(x_+,\pi _+,x_a,\pi _a)`$.

Now for a system with $`n`$ degrees of freedom and $`r`$ first class constraints $`\varphi ^\alpha `$, the matrix element of the $`S`$ \- matrix is given by Faddeev and Popov  as

$$Out|S|In=\underset{t}{}d\mu (q_j,p_j)\mathrm{exp}[i\{_{\mathrm{}}^{\mathrm{}}𝑑t(p_j\dot{q_j}H_0)\}],j=1,\mathrm{},n,$$
(44)
where the measure of integration is given as

$$d\mu =det|\{\varphi ^\alpha ,\chi ^\beta \}|\underset{\alpha =1}{\overset{r}{}}\delta (\chi ^\alpha )\delta (\varphi ^\alpha )\underset{j=1}{\overset{n}{}}dq^jdp_j,$$
(45)
and $`\chi _a`$ are $`r`$\- gauge constraints.

If we perform now the usual path integral quantization  using (44) for the system (13), one must choose one time dependent gauge fixing (equation (25)) condition to obtain the path integral quantization over the canonical phase-space coordinates.

## 4  Conclusion

Path integral quantization for relativistic particle in plane wave is obtained by using the canonical path integral formulation \[8-10\]. In this approach, since the integrability conditions $`dH_{}^{}{}_{\tau }{}^{}=0`$ and $`dH_{}^{}{}_{x_{}}{}^{}=0`$ are satisfied identically, this system is integrable. Hence, the canonical phase-space coordinates $`(x_+,\pi _+,x_a,\pi _a)`$ are obtained in terms of parameter $`(x_{})`$. In this case the path integral, then follows directly as given in (43) without using any gauge fixing conditions.

When applying the analog of Faddeev’s method  to this model, one may choose gauge fixing of the form $`x_{}\tau =0`$, so that one can integrate over the extended phase space coordinates $`(x_{},\pi _{},x_+,\pi _+,x_a,\pi _a)`$ and after integration over the redundant variables $`(x_{},\pi _{})`$, one can arrive at the result (43). Besides different choice of gauge fixing will lead to different quantum theories. In this case there are many measures one can use to define an integral. In other words considerable care must be taken with the choice of gauge fixing conditions in order to arrive at the correct result.

As a conclusion, the path integral for relativistic particle in a plane wave is obtained without using any gauge fixing conditions and without the problems of the measure. Besides the canonical action integral is obtained with the consequncent constraints (30,31) are enforced explicitly without Lagrange multipliers.


format: text

text:
# Fermion helicity flip by parity violating torsion

## Acknowledgements

We acknowledge stimulating discussion with A. Chatterjee, T. De, R. Gandhi, B. Mukhopadhyay and A. Raychaudury. S.S. is supported by project grant no. 98/37/16/BRNS cell/676 from the Board Of Research In Nuclear Sciences, Department Of Atomic Energy, Government Of India. A.S. is supported by a scholarship from Cambridge Commonwealth Trust.


text:
# Implementation of analytical Hartree-Fock gradients for periodic systems

## I Introduction

The determination of equilibrium structure is one of the most important targets in electronic structure calculations. In surface science especially, theoretical calculations of surface structures are of high importance to explain and support experimental results. Therefore, a fast structural optimization is an important issue in modern electronic structure codes. Finding minima in energy surfaces is substantially simplified by the availability of analytical gradients. As a rule of thumb, availability of analytical gradients improves the efficiency by a factor of order $`N`$ with $`N`$ being the number of parameters to be optimized. UK’s Collaborative Computational Project 3 has therefore supported the implementation of analytical gradients in the electronic structure code CRYSTAL . This implementation will also be valuable for future projects which require analytical gradients as a prerequisite. Another advantage of having analytical gradients is that higher derivatives can be obtained with less numerical noise (e.g. the 2nd derivative has less numerical noise when only one numerical differentiation is necessary).

CRYSTAL is capable of performing Hartree-Fock and density-functional calculations for systems with any periodicity (i.e. molecules, polymers, slabs and solids). The periodicity is ”cleanly” implemented in the sense that, for example, a slab is considered as an object periodic in two dimensions and is not repeated in the third dimension with one slab being separated from the others by vacuum layers. The code is based on Gaussian type orbitals and the technology is therefore in many parts similar to that of molecular quantum chemistry codes. As the density-functional part of the code relies in big parts on the Hartree-Fock part, the strategy of the project was to implement Hartree-Fock gradients first.

The implementation of Hartree-Fock gradients for multicenter basis sets was pioneered by Pulay; the theory had already been derived earlier independently . Meanwhile, analytical gradients have been implemented in many molecular codes, and several review articles have appeared (see, e.g., references ).

Substantial work has also been done in the case of one-dimensional periodicity: Hartree-Fock gradients with respect to nuclear coordinates and with respect to the lattice vector have already been implemented in codes periodic in one dimension . Moreover, correlated calculations based on the MP2 scheme  and MP2 gradients  have been coded. Also, density functional gradients have been implemented . Even second derivatives at the Hartree-Fock level have meanwhile been coded .

The implementation of Hartree-Fock gradients with respect to nuclear coordinates in CRYSTAL is to the best of our knowledge the first implementation for the case of 2- and 3-dimensional periodicity. The aim of this article is to describe the implementation of the gradients in the code, with an emphasis on the technical aspects. Therefore, the article is supposed to complement our first article on the purely theoretical aspects . An attempt of a detailed description is made; however, as the whole code is undergoing constant changes, it can not be too detailed. For example, it did not seem advisable to give any variable names because they have already undergone major changes after the code moved to Fortran 90 with the possibility of longer variable names.

The article is structured as follows: In section II, we give a brief introduction to Gaussian and Hermite Gaussian type basis functions. The definition of the density matrix is given in section III. The individual integrals, their derivatives, and details of the implementation are discussed in section IV. Formulas for total energy and gradient are given in section V. The structure of the gradient code is explained in section VI, followed by examples in section VII and the conclusion.

## II Basis functions

Two sets of basis functions are relevant for CRYSTAL: firstly, unnormalized spherical Gaussian type functions, in a polar coordinate system characterized by the set of variables $`(|\stackrel{}{r}|,\vartheta ,\phi )`$, and centered at $`\stackrel{}{A}`$. They are defined as

$$S(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,m)=|\stackrel{}{r}\stackrel{}{A}|^{2n+l}\mathrm{P}_l^{|m|}(\mathrm{cos}\vartheta )\mathrm{exp}(\mathrm{i}m\phi )\mathrm{exp}(\alpha |\stackrel{}{r}\stackrel{}{A}|^2)$$
(1)
with $`\mathrm{P}_l^{|m|}`$ being the associated Legendre function. CRYSTAL uses real spherical Gaussian type functions defined as

$`R(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,0)=S(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,0)`$
$`R(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,|m|)=\mathrm{Re}S(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,|m|)`$
$`R(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,|m|)=\mathrm{Im}S(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,|m|)`$
This is in the following denoted as $`\varphi _\mu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\mu ,n,l,m)=N_\mu R(\alpha ,\stackrel{}{r}\stackrel{}{A}_\mu ,n,l,m)`$, with the normalization $`N_\mu `$. $`\mu `$ is an index enumerating the basis functions in the reference cell (e.g. the primitive unit cell). In fact, CRYSTAL uses only basis functions with quantum number $`n=0`$ and angular momentum $`l`$=0,1 or 2 (i.e. $`s`$, $`p`$ or $`d`$ functions).

The exponents are defined by the user of the code. A huge amount of basis sets for molecular calculations is available in the literature and on the world wide web; also for periodic systems a large number of basis sets has been published. Molecular basis sets can, with a little effort, be adopted for solid state calculations. High exponents which are used to describe core electrons do not have to be adjusted, but exponents with low values (e.g. less than 1 $`a_0^2`$, with $`a_0`$ being the Bohr radius) should be reoptimized for the corresponding solid. Very diffuse exponents should be omitted because they cause linear dependence problems in periodic systems.

A second type of basis functions, which CRYSTAL uses internally to evaluate the integrals, is the Hermite Gaussian type function (HGTF) which is defined as:

$`\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{A},t,u,v)=\left({\displaystyle \frac{}{A_x}}\right)^t\left({\displaystyle \frac{}{A_y}}\right)^u\left({\displaystyle \frac{}{A_z}}\right)^v\mathrm{exp}(\gamma |\stackrel{}{r}\stackrel{}{A}|^2)`$ (2)
CRYSTAL uses the McMurchie-Davidson algorithm to evaluate the integrals. The basic idea of this algorithm is to map the product of two spherical Gaussian type functions on two centers onto a set of Hermite Gaussian type functions at one center.

$`S(\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{B},\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})S(\alpha ,\stackrel{}{r}\stackrel{}{A},n,l,m)=`$ (3)
$`{\displaystyle \underset{t,u,v}{}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)`$ (4)
with $`\gamma =\alpha +\stackrel{~}{\alpha }`$ and $`\stackrel{}{P}=\frac{\alpha \stackrel{}{A}+\stackrel{~}{\alpha }\stackrel{}{B}}{\alpha +\stackrel{~}{\alpha }}`$.

The starting point $`E(0,0,0,0,0,0,0,0,0)=\mathrm{exp}(\frac{\alpha \stackrel{~}{\alpha }}{\alpha +\stackrel{~}{\alpha }}|\stackrel{}{B}\stackrel{}{A}|^2)`$ is derived from the Gaussian product rule :

$`\mathrm{exp}(\alpha |\stackrel{}{r}\stackrel{}{A}|^2)\mathrm{exp}(\stackrel{~}{\alpha }|\stackrel{}{r}\stackrel{}{B}|^2)=\mathrm{exp}\left({\displaystyle \frac{\alpha \stackrel{~}{\alpha }}{\alpha +\stackrel{~}{\alpha }}}|\stackrel{}{B}\stackrel{}{A}|^2\right)\mathrm{exp}\left((\alpha +\stackrel{~}{\alpha })\left|\stackrel{}{r}{\displaystyle \frac{\alpha \stackrel{}{A}+\stackrel{~}{\alpha }\stackrel{}{B}}{\alpha +\stackrel{~}{\alpha }}}\right|^2\right)`$ (5)
As indicated in section IV, all the integrals can be expressed with the help of the coefficients $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)`$ . These coefficients are generated by recursion relations . They are zero for the case $`t+u+v>2n+2\stackrel{~}{n}+l+\stackrel{~}{l}`$ and for all negative values of $`t,u`$ or $`v`$. CRYSTAL uses only basis functions with $`n=0`$. Therefore, there are $`\frac{(l+\stackrel{~}{l}+1)(l+\stackrel{~}{l}+2)(l+\stackrel{~}{l}+3)}{3!}`$ coefficients $`E(0,\stackrel{~}{l},\stackrel{~}{m},0,l,m,t,u,v)`$ for fixed values of $`l,m,\stackrel{~}{l},\stackrel{~}{m}`$. As the maximum angular quantum number is $`l=2`$, this results in 25 possible combinations of $`m`$ and $`\stackrel{~}{m}`$. Therefore, the maximum number of coefficients is $`25\times 35=875`$. These coefficients are pre-programmed in the subroutine DFAC3. Pre-programming is the fastest possible way of evaluating these coefficients which is important because this is one of the key issues of the integral calculation. On the other hand, the code has become inflexible as no $`E`$-coefficients are available for higher quantum numbers.

Derivatives of Gaussian type functions are again Gaussian type functions. Therefore, the evaluation of gradients is closely related to the evaluation of integrals. In a similar way as all the integrals can be expressed with the help of coefficients $`E`$, all the derivatives of the integrals can be expressed with the help of coefficients for the gradients, $`G_x^A,G_y^A,G_z^A`$. These $`G`$-coefficients can be obtained with recursion relations derived by Saunders . The recursions are similar to the ones for the $`E`$-coefficients. However, as the existing subroutine DFAC3 cannot compute the $`G`$-coefficients, the recursions were newly coded. This has in addition the advantage that, by small modifications of the new subroutines, $`E`$-coefficients for higher quantum numbers than $`l=\stackrel{~}{l}=2`$ can now be computed by recursion. There are three sets of $`G`$-coefficients because of the three spatial directions. The $`G`$-coefficients are zero for the case $`t+u+v>2n+2\stackrel{~}{n}+l+\stackrel{~}{l}+1`$ and for all negative values of $`t,u`$ or $`v`$. This means that for a maximum quantum number of $`l=2`$, there are $`3\times 5\times 5\times 56=4200`$ coefficients. Three other sets of $`G`$-coefficients are necessary because of the second center. However, the sets on the second center are closely related to the sets on the first center and can be derived from them in an efficient way .

## III Density matrix

After solving the Hartree-Fock equations , the crystalline orbitals are linear combinations of Bloch functions

$$\mathrm{\Psi }_i(\stackrel{}{r},\stackrel{}{k})=\underset{\mu }{}a_{\mu i}(\stackrel{}{k})\psi _\mu (\stackrel{}{r},\stackrel{}{k})$$
(6)
which are expanded in terms of real spherical Gaussian type functions

$$\psi _\mu (\stackrel{}{r},\stackrel{}{k})=N_\mu \underset{\stackrel{}{g}}{}R(\alpha ,\stackrel{}{r}\stackrel{}{A}_\mu \stackrel{}{g},n,l,m)\mathrm{e}^{\mathrm{i}\stackrel{}{k}\stackrel{}{g}}$$
(7)
The sum over $`\stackrel{}{g}`$ is over all direct lattice vectors.

In the case of closed shell, spin-restricted Hartree-Fock, the spin-free density matrix in reciprocal space is defined as

$`P_{\mu \nu }(\stackrel{}{k})=2{\displaystyle \underset{i}{}}a_{\mu i}(\stackrel{}{k})a_{\nu i}^{}(\stackrel{}{k})\mathrm{\Theta }(ϵ_Fϵ_i(\stackrel{}{k}))`$ (8)
with the Fermi energy $`ϵ_F`$ and the Heaviside function $`\mathrm{\Theta }`$; $`i`$ is an index enumerating the eigenvalues.

In the case of unrestricted Hartree-Fock (UHF) , we use the notation

$$\mathrm{\Psi }_i^{}(\stackrel{}{r},\stackrel{}{k})=\underset{\mu }{}a_{\mu i}^{}(\stackrel{}{k})\psi _\mu (\stackrel{}{r},\stackrel{}{k})$$
(9)
and

$$\mathrm{\Psi }_i^{}(\stackrel{}{r},\stackrel{}{k})=\underset{\mu }{}a_{\mu i}^{}(\stackrel{}{k})\psi _\mu (\stackrel{}{r},\stackrel{}{k})$$
(10)
for the crystalline orbitals with up and down spin, respectively. We define the density matrices

$`P_{\mu \nu }^{}(\stackrel{}{k})={\displaystyle \underset{i}{}}a_{\mu i}^{}(\stackrel{}{k})a_{\nu i}^{}(\stackrel{}{k})\mathrm{\Theta }(ϵ_Fϵ_i^{}(\stackrel{}{k}))`$ (11)
for up spin and

$`P_{\mu \nu }^{}(\stackrel{}{k})={\displaystyle \underset{i}{}}a_{\mu i}^{}(\stackrel{}{k})a_{\nu i}^{}(\stackrel{}{k})\mathrm{\Theta }(ϵ_Fϵ_i^{}(\stackrel{}{k}))`$ (12)
for down spin. In the following, $`P_{\mu \nu }`$ refers to the sum $`P_{\mu \nu }^{}+P_{\mu \nu }^{}`$ in the UHF case.

The density matrices in real space $`P_{\mu \stackrel{}{0}\nu \stackrel{}{g}},P_{\mu \stackrel{}{0}\nu \stackrel{}{g}}^{},P_{\mu \stackrel{}{0}\nu \stackrel{}{g}}^{}`$ are obtained by Fourier transformation.

## IV Integrals and their derivatives

The calculation of the integrals is fundamental to all quantum chemistry programs. CRYSTAL uses two integral packages: a package derived from GAUSSIAN70  is the default for calculations when only $`s`$ and $`sp`$ shells are used; alternatively Saunders’ ATMOL Gaussian integral package can be used and it must be used for cases when $`p`$ or $`d`$ functions are involved. The implementation of gradients has been done with routines based on the ATMOL package. This is not a restriction, and it is possible to use routines based on GAUSSIAN70 for the integrals and routines based on ATMOL for the gradients.

The calculation of the integrals is essentially controlled from MONMAD and MONIRR for one-electron integrals and from SHELLC or SHELLX for the bielectronic integrals. SHELLC is used in the case of non-direct SCF, i.e. when the integrals are written to disk and read in each cycle. SHELLX is the direct version when the integrals are computed in each cycle without storing them on disk. The direct mode is the preferred one when the integral file is too big or when input/output to disk is too slow. The gradients are computed only once after the last iteration, when convergence is achieved. Therefore, a direct implementation of gradients has been done.

One of the bottlenecks of the CRYSTAL code is the restriction to a highest quantum number of $`l=2`$, i.e. the code can only cope with $`s`$, $`p`$, $`sp`$ and $`d`$ functions, but not with basis functions with higher angular momentum. Introducing gradients, however, is similar to increasing the quantum number from $`d`$ to $`f`$ for the corresponding basis function. This means that many subroutines had to be extended to higher quantum numbers, and array dimensions in the whole code had to be adjusted.

### A One-electron integrals

In this section we summarize the appearing types of integrals and the corresponding gradients. We restrict the description to the x-component of the gradient; y- and z-component can be obtained in similar way. Note that the integrals depend on the dimension because of the Ewald scheme used. Therefore, there are four different routines for the one-electron integrals for the case of 0,1,2 and 3-dimensional periodicity: CJAT0, CJAT1, CJAT2 and CJAT3. Similarly, four gradient routines have been developed which have been given the preliminary names CJAT0G, CJAT1G, CJAT2G and CJAT3G. These routines calculate all the one-electron integrals except for the multipolar integrals which are computed in POLIPA (with the corresponding gradient routine OLIPAG).

#### 1 Overlap integral

The basic integral is the overlap integral:

$`S_{\mu \stackrel{}{g_1}\nu \stackrel{}{g_2}}={\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu \stackrel{}{g}_1,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g}_2,n,l,m)\mathrm{d}^3\mathrm{r}}=`$ (13)
$`{\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}=`$ (14)
$`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,0)\left({\displaystyle \frac{\pi }{\gamma }}\right)^{\frac{3}{2}}`$ (15)
The x-component of the gradient with respect to center $`A_\mu `$ is obtained as

$`{\displaystyle \frac{}{A_{\mu ,x}}}S_{\mu \stackrel{}{g_1}\nu \stackrel{}{g_2}}=`$ (16)
$`{\displaystyle \frac{}{A_{\mu ,x}}}{\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu \stackrel{}{g}_1,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g}_2,n,l,m)\mathrm{d}^3\mathrm{r}}=`$ (17)
$`{\displaystyle \frac{}{A_{\mu ,x}}}{\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}=`$ (18)
$`{\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}=`$ (19)
$`G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,0)\left({\displaystyle \frac{\pi }{\gamma }}\right)^{\frac{3}{2}}`$ (20)
Equation 20 thus defines the coefficients $`G_x^{A_\mu }`$; similarly the coefficients $`G_y^{A_\mu },G_z^{A_\mu },G_x^{A_\nu },G_y^{A_\nu },G_z^{A_\nu }`$ can be defined.

In the following, we use the identity

$`S_{\mu \stackrel{}{g_1}\nu \stackrel{}{g_2}}=S_{\mu \stackrel{}{0}\nu (\stackrel{}{g_2}\stackrel{}{g_1})}=S_{\mu \stackrel{}{0}\nu \stackrel{}{g}}`$.

#### 2 Kinetic energy integrals

In equation 21, the expression for the kinetic energy integrals for the case of spherical Gaussian type functions is reiterated :

$`T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}=`$ (21)
$`{\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\left(\frac{1}{2}\mathrm{\Delta }_\stackrel{}{r}\right)\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)\mathrm{d}^3\mathrm{r}}=`$ (22)
$`n(2n+2l+1){\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n1,l,m)\mathrm{d}^3\mathrm{r}}+`$ (23)
$`\alpha (4n+2l+3){\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)\mathrm{d}^3\mathrm{r}}`$ (24)
$`2\alpha ^2{\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n+1,l,m)\mathrm{d}^3\mathrm{r}}=`$ (25)
$`n(2n+2l+1){\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n1,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}+`$ (26)
$`\alpha (4n+2l+3){\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}`$ (27)
$`2\alpha ^2{\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n+1,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}`$ (28)
The x-component of the gradient is therefore:

$`{\displaystyle \frac{}{A_{\mu ,x}}}T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}=`$ (29)
$`n(2n+2l+1){\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n1,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}+`$ (30)
$`\alpha (4n+2l+3){\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}`$ (31)
$`2\alpha ^2{\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n+1,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\mathrm{d}^3\mathrm{r}}`$ (32)
As CRYSTAL uses spherical Gaussian type functions with $`n=0`$, this reduces to

$`{\displaystyle \frac{}{A_{\mu ,x}}}T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}=`$ (33)
$`\left({\displaystyle \frac{\pi }{\gamma }}\right)^{\frac{3}{2}}\alpha (4n+2l+3)G_x^{A_\mu }(0,\stackrel{~}{l},\stackrel{~}{m},0,l,m,0,0,0)`$ (34)
$`2\left({\displaystyle \frac{\pi }{\gamma }}\right)^{\frac{3}{2}}\alpha ^2G_x^{A_\mu }(0,\stackrel{~}{l},\stackrel{~}{m},1,l,m,0,0,0)`$ (35)
Explicit differentiation with respect to the other center $`\stackrel{}{A}_\nu `$ is more difficult because the kinetic energy operator applies to that center. However, the differentiation can easily be avoided by applying translational invariance:

$`{\displaystyle \frac{}{A_{\mu ,x}}}T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}={\displaystyle \frac{}{A_{\nu ,x}}}T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}`$ (36)
#### 3 Nuclear attraction integrals

The nuclear attraction integrals are defined as

$`N_{\mu \stackrel{}{0}\nu \stackrel{}{g}}={\displaystyle \underset{a}{}}Z_a{\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})A(\stackrel{}{r}\stackrel{}{A}_a)\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)\mathrm{d}^3\mathrm{r}}=`$ (37)
$`{\displaystyle \underset{a}{}}Z_a{\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)A(\stackrel{}{r}\stackrel{}{A}_a)\mathrm{d}^3\mathrm{r}}`$ (38)
where $`A`$ is the Coulomb potential function in the molecular case, the Euler-MacLaurin potential function for systems periodic in one dimension , Parry’s potential function  for systems periodic in two dimensions, and the Ewald potential function for systems periodic in three dimensions . The summation with respect to $`a`$ runs over all nuclei of the primitive unit cell.

The x-component of the partial derivative with respect to the center $`A_{\mu ,x}`$ is obtained as:

$`{\displaystyle \frac{}{A_{\mu ,x}}}N_{\mu \stackrel{}{0}\nu \stackrel{}{g}}=`$ (39)
$`{\displaystyle \underset{a}{}}Z_a{\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)A(\stackrel{}{r}\stackrel{}{A}_a)\mathrm{d}^3\mathrm{r}}`$ (40)
In the same way, the partial derivative with respect to $`A_{\nu ,x}`$ is obtained. The partial derivative with respect to the set of third centers $`\stackrel{}{A}_a`$ is obtained by translational invariance: for each center $`\stackrel{}{A}_a`$, there is a derivative with value
$`\frac{}{\stackrel{}{A}_\mu }\frac{}{\stackrel{}{A}_\nu }`$.

#### 4 Multipolar integrals

The electronic charge density is expressed with a lattice basis as:

$`\rho (\stackrel{}{r})={\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}\varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)`$ (41)
Then, the Ewald potential due to this charge density is given by:

$`\mathrm{\Phi }^{ew}(\rho ;\stackrel{}{r})={\displaystyle A(\stackrel{}{r}\stackrel{}{r}^{})\rho (\stackrel{}{r}^{})d^3r^{}}`$ (42)
The Ewald energy of the electons (i.e. the Ewald energy of the electrons in the primitive unit cell with all the electrons) is obtained as

$`E={\displaystyle \frac{1}{2}}{\displaystyle \rho (\stackrel{}{r})A(\stackrel{}{r}\stackrel{}{r}^{})\rho (\stackrel{}{r}^{})d^3rd^3r^{}}`$ (43)
For efficiency reasons, the calculation of the Ewald potential is done approximatively. A multipolar expansion up to an order $`L`$ is performed for the charge distribution in the long range. Therefore, the electrons do not feel the Ewald potential created by the correct charge distribution, but the Ewald potential created by the multipolar moments. It is thus necessary to compute the multipolar moments of the charge distribution which are defined as

$`\eta _l^m(\rho _c;\stackrel{}{A}_c)={\displaystyle \rho _c(\stackrel{}{r})X_l^m(\stackrel{}{r}\stackrel{}{A}_c)\mathrm{d}^3\mathrm{r}}`$ (44)
with $`X_l^m`$ being regular solid harmonics  and the charge $`\rho _c(\stackrel{}{r})`$ defined as

$`\rho _c(\stackrel{}{r})={\displaystyle \underset{\stackrel{}{g},\mu c,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}\varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)=`$ (45)
$`{\displaystyle \underset{\stackrel{}{g},\mu c,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \underset{t,u,v}{}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)`$ (46)
$`c`$ is an index for the shell. The total electronic charge $`\rho (\stackrel{}{r})`$ is thus obtained by summing over all shells $`c`$:

$`\rho (\stackrel{}{r})={\displaystyle \underset{c}{}}\rho _c(\stackrel{}{r})`$ (47)
In CRYSTAL, the multipole is located at center $`\stackrel{}{A}_\mu `$ and therefore it is convenient to take the derivative with respect to center $`\stackrel{}{A}_\nu `$ and obtain the derivative with respect to $`\stackrel{}{A}_\mu `$ by translational invariance. The expression computed for the gradients is thus

$`{\displaystyle \underset{\stackrel{}{g},\mu c,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \frac{}{A_{\nu ,x}}\left(\varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)X_l^m(\stackrel{}{r}\stackrel{}{A}_\mu )\right)\mathrm{d}^3\mathrm{r}}=`$ (48)
$`{\displaystyle \underset{\stackrel{}{g},\mu c,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \underset{t,u,v}{}G_x^{A_\nu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)X_l^m(\stackrel{}{r}\stackrel{}{A}_\mu )\mathrm{d}^3\mathrm{r}}`$ (49)
#### 5 Field integrals

If the electronic charge distribution is approximated with an expansion up to the maximum quantum number $`L`$, the Ewald potential of this model charge distribution is obtained as

$`\mathrm{\Phi }^{ew}(\rho ^{model};\stackrel{}{r})={\displaystyle \underset{c}{}}\mathrm{\Phi }^{ew}(\rho _c^{model};\stackrel{}{r})={\displaystyle \underset{c}{}}{\displaystyle \underset{l=0}{\overset{L}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}\eta _l^m(\rho _c;\stackrel{}{A}_c)Z_l^m(\widehat{\stackrel{}{A}_c})A(\stackrel{}{r}\stackrel{}{A}_c)`$ (50)
with $`Z_l^m(\widehat{\stackrel{}{A}_c})`$ being the spherical gradient operator in a renormalized form. The model charge distribution is expressed as

$`\rho _c^{\mathrm{model}}(\stackrel{}{r})={\displaystyle \underset{l=0}{\overset{L}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}\eta _l^m(\rho _c;\stackrel{}{A}_c)\delta _l^m(\stackrel{}{A}_c,\stackrel{}{r})`$ (51)
and

$`\delta _l^m(\stackrel{}{A_c},\stackrel{}{r})=\underset{\alpha \mathrm{}}{lim}Z_l^m(\widehat{\stackrel{}{A}_c})\mathrm{\Lambda }(\alpha ,\stackrel{}{r}\stackrel{}{A}_c,0,0,0)`$ (52)
The integral of the electronic charge distribution and the Ewald potential function is required which gives rise to the field integrals which are defined as follows:

$`M_{l\mu \stackrel{}{0}\nu \stackrel{}{g}c}^m=`$ (53)
$`Z_l^m(\widehat{\stackrel{}{A}_c}){\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)\left[A(\stackrel{}{r}\stackrel{}{A}_c)\underset{\stackrel{}{n}}{\overset{pen}{}}\frac{1}{|\stackrel{}{r}\stackrel{}{A}_c\stackrel{}{n}|}\right]\mathrm{d}^3\mathrm{r}}=`$ (54)
$`Z_l^m(\widehat{\stackrel{}{A}_c}){\displaystyle \underset{t,u,v}{}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\left[A(\stackrel{}{r}\stackrel{}{A}_c)\underset{\stackrel{}{n}}{\overset{pen}{}}\frac{1}{|\stackrel{}{r}\stackrel{}{A}_c\stackrel{}{n}|}\right]\mathrm{d}^3\mathrm{r}}`$ (55)
The term $`\left[A(\stackrel{}{r}\stackrel{}{A}_c)_\stackrel{}{n}^{pen}\frac{1}{|\stackrel{}{r}\stackrel{}{A}_c\stackrel{}{n}|}\right]`$ instead of $`A(\stackrel{}{r}\stackrel{}{A}_c)`$ appears because the multipolar approximation is only done for the charge distribution in the long range. The penetration depth $`pen`$ is a certain threshold for which the integrals are evaluated exactly .

For the gradients, the derivative with respect to all the centers is needed. The partial derivative with respect to $`A_{\mu ,x}`$ is obtained as

$`{\displaystyle \frac{}{A_{\mu ,x}}}M_{l\mu \stackrel{}{0}\nu \stackrel{}{g}c}^m=`$ (56)
$`Z_l^m(\widehat{\stackrel{}{A}_c}){\displaystyle \underset{t,u,v}{}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\left[A(\stackrel{}{r}\stackrel{}{A}_c)\underset{\stackrel{}{n}}{\overset{pen}{}}\frac{1}{|\stackrel{}{r}\stackrel{}{A}_c\stackrel{}{n}|}\right]\mathrm{d}^3\mathrm{r}}`$ (57)
In similar way, the partial derivative with respect to center $`\stackrel{}{A}_\nu `$ is computed. Finally, the partial derivatives with respect to the centers $`\stackrel{}{A}_c`$ are obtained from translational invariance.

#### 6 Spheropole

This term arises because the charge distribution is approximated by a model charge distribution in the long range:

$`\mathrm{\Phi }^{ew}(\rho _c;\stackrel{}{r})=\mathrm{\Phi }^{ew}(\rho _c^{model};\stackrel{}{r})+\mathrm{\Phi }^{ew}(\rho _c\rho _c^{model};\stackrel{}{r})=\mathrm{\Phi }^{ew}(\rho _c^{model};\stackrel{}{r})+\mathrm{\Phi }^{coul}(\rho _c\rho _c^{model};\stackrel{}{r})+Q_c`$ (58)
The calculation of the Coulomb potential $`\mathrm{\Phi }^{coul}(\rho _c\rho _c^{model};\stackrel{}{r})`$ is restricted to contributions from those charges inside the penetration depth $`pen`$. The use of the Coulomb potential $`\mathrm{\Phi }^{coul}(\rho _c\rho _c^{model};\stackrel{}{r})`$ instead of the Ewald potential $`\mathrm{\Phi }^{ew}(\rho _c\rho _c^{model};\stackrel{}{r})`$ is correct, if $`\rho _c\rho _c^{model}`$ is of zero charge, dipole, quadrupole and spherical second moment. However, this condition leads to a correction in the three-dimensional case: although the difference $`\rho _c\rho _c^{model}`$ has zero charge, dipole and quadrupole moment, it has in general a non-zero spherical second moment $`Q_c`$. Therefore, the potential must be shifted by $`Q`$ defined as:

$`Q={\displaystyle \underset{c}{}}Q_c={\displaystyle \underset{c}{}}{\displaystyle \frac{2\pi }{3V}}{\displaystyle (\rho _c(\stackrel{}{r})\rho _c^{\mathrm{model}}(\stackrel{}{r}))|\stackrel{}{r}|^2\mathrm{d}^3\mathrm{r}}`$ (59)
Three types of contributions are obtained : zero, first and second order HGTFs. They have to be combined with the corresponding $`E`$-coefficient. For the zeroth order, a contribution of

$`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,0)\left(\frac{3}{2\gamma }+\left(\stackrel{}{A}_\mu \stackrel{}{P}\right)^2\right)`$

is computed. The derivative is therefore

$`{\displaystyle \frac{}{A_{\mu ,x}}}\left(E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,0)\left({\displaystyle \frac{3}{2\gamma }}+(\stackrel{}{A}_\mu \stackrel{}{P})^2\right)\right)`$ (60)
To obtain the derivative $`\frac{}{A_{\mu ,x}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,0)`$, we use the identity

$`{\displaystyle \frac{}{A_{\mu ,x}}}\left({\displaystyle \underset{t,u,v}{}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)\right)=`$ (61)
$`{\displaystyle \underset{t,u,v}{}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v){\displaystyle \frac{\stackrel{~}{\alpha }}{\gamma }}\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t+1,u,v)+`$ (62)
$`{\displaystyle \underset{t,u,v}{}}\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v){\displaystyle \frac{}{A_{\mu ,x}}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)=`$ (63)
$`{\displaystyle \underset{t,u,v}{}}G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v)`$ (64)
which gives

$`{\displaystyle \frac{}{A_{\mu ,x}}}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v)=G_x^{A_\mu }(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t,u,v){\displaystyle \frac{\stackrel{~}{\alpha }}{\gamma }}E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,t1,u,v)`$ (65)
A similar operation is necessary for the components with $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,1,0,0)`$, $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,1,0)`$ and $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,1)`$ (first order HGTFs) which are multiplied with prefactors $`2(P_xA_{\mu ,x})`$, $`2(P_yA_{\mu ,y})`$ and $`2(P_zA_{\mu ,z})`$, respectively. Finally, derivatives of the products of $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,2,0,0)`$, $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,2,0)`$ and $`E(\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m},n,l,m,0,0,2)`$ (second order HGTFs) with 2 are required.

### B Bielectronic integrals

We define a bielectronic integral as

$`B_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{n}\sigma \stackrel{}{n}+\stackrel{}{h}}=`$ (66)
$`{\displaystyle \frac{\varphi _\mu (\alpha _1,\stackrel{}{r}\stackrel{}{A}_\mu ,n_1,l_1,m_1)\varphi _\nu (\alpha _2,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n_2,l_2,m_2)}{|\stackrel{}{r}\stackrel{}{r}^{}|}}`$ (67)
$`\varphi _\tau (\alpha _3,\stackrel{}{r}^{}\stackrel{}{A}_\tau \stackrel{}{n},n_3,l_3,m_3)\varphi _\sigma (\alpha _4,\stackrel{}{r}^{}\stackrel{}{A}_\sigma \stackrel{}{n}\stackrel{}{h},n_4,l_4,m_4)\mathrm{d}^3\mathrm{r}\mathrm{d}^3\mathrm{r}^{}=`$ (68)
$`{\displaystyle \underset{t,u,v}{}}E(n_1,l_1,m_1,n_2,l_2,m_2,t,u,v){\displaystyle \underset{t^{},u^{},v^{}}{}}E(n_3,l_3,m_3,n_4,l_4,m_4,t^{},u^{},v^{})[t,u,v|{\displaystyle \frac{1}{|\stackrel{}{r}\stackrel{}{r}^{}|}}|t^{},u^{},v^{}]`$ (69)
The expression $`[t,u,v|\frac{1}{|\stackrel{}{r}\stackrel{}{r}^{}|}|t^{},u^{},v^{}]`$ is defined as

$`[t,u,v|{\displaystyle \frac{1}{|\stackrel{}{r}\stackrel{}{r}^{}|}}|t^{},u^{},v^{}]=`$ (70)
$`{\displaystyle }{\displaystyle }\mathrm{\Lambda }(\gamma ,\stackrel{}{r}\stackrel{}{P},t,u,v){\displaystyle \frac{1}{|\stackrel{}{r}\stackrel{}{r}^{}|}}\mathrm{\Lambda }(\gamma ^{},\stackrel{}{r}^{}\stackrel{}{P}^{},t,^{}u^{},v^{})\mathrm{d}^3\mathrm{r}\mathrm{d}^3\mathrm{r}^{}`$ (71)
The partial derivative with respect to $`A_{\mu ,x}`$ is obtained as

$`{\displaystyle \frac{}{A_{\mu ,x}}}B_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{n}\sigma \stackrel{}{n}+\stackrel{}{h}}=`$ (72)
$`{\displaystyle \underset{t,u,v}{}}G_x^{A_\mu }(n_1,l_1,m_1,n_2,l_2,m_2,t,u,v){\displaystyle \underset{t^{},u^{},v^{}}{}}E(n_3,l_3,m_3,n_4,l_4,m_4,t^{},u^{},v^{})[t,u,v|{\displaystyle \frac{1}{|\stackrel{}{r}\stackrel{}{r}^{}|}}|t^{},u^{},v^{}]`$ (73)
Similarly, gradients with respect to the other centers are obtained. One of the gradients can be obtained by translational invariance if the other three gradients have been computed.

In the context of periodic systems, it is necessary to perform summations over the lattice vectors $`\stackrel{}{g},\stackrel{}{h},\stackrel{}{n}`$. We define a Coulomb integral as follows

$`C_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}={\displaystyle \underset{\stackrel{}{n}}{\overset{pen}{}}}B_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{n}\sigma \stackrel{}{n}+\stackrel{}{h}}`$ (74)
Similarly, we define an exchange integral as follows:

$`X_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}={\displaystyle \underset{\stackrel{}{n}}{}}B_{\mu \stackrel{}{0}\tau \stackrel{}{n}\nu \stackrel{}{g}\sigma \stackrel{}{n}+\stackrel{}{h}}`$ (75)
## V Total energy and gradient

### A Total energy

The correct summation of the Coulomb energy is the most severe problem of the total energy calculation. The individual contributions to the Coulomb energy, such as for example the nuclear-nuclear interaction, are divergent for periodic systems. Thus, a scheme based on the Ewald method is used to sum the individual contributions . The total energy is then expressed as the sum of kinetic energy $`E^{kin}`$, the Ewald energies of the nuclear-nuclear repulsion $`E^{NN}`$, nuclear-electron attraction $`E^{coulnuc}`$, electron-electron repulsion $`E^{coulel}`$, and finally the exchange energy $`E^{exchel}`$.

$`E^{\mathrm{total}}=E^{\mathrm{kinetic}}+E^{\mathrm{NN}}+E^{\mathrm{coul}\mathrm{nuc}}+E^{\mathrm{coul}\mathrm{el}}+E^{\mathrm{exch}\mathrm{el}}=`$ (76)
$`={\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}+E^{\mathrm{NN}}`$ (77)
$`{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \underset{a}{}}Z_a{\displaystyle \varphi _\mu (\stackrel{~}{\alpha },\stackrel{}{r}\stackrel{}{A}_\mu ,\stackrel{~}{n},\stackrel{~}{l},\stackrel{~}{m})\varphi _\nu (\alpha ,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n,l,m)A(\stackrel{}{r}\stackrel{}{A}_a)\mathrm{d}^3\mathrm{r}}`$ (78)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}\left(QS_{\mu \stackrel{}{0}\nu \stackrel{}{g}}+{\displaystyle \underset{\stackrel{}{h},\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}C_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}{\displaystyle \underset{c}{}}{\displaystyle \underset{l=0}{\overset{L}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}\eta _l^m(\rho _c;\stackrel{}{A}_c)M_{l\mu \stackrel{}{0}\nu \stackrel{}{g}c}^m\right)`$ (79)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}^{}{\displaystyle \underset{\stackrel{}{h},\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}^{}X_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}^{}{\displaystyle \underset{\stackrel{}{h},\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}^{}X_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}`$ (80)
### B Gradient of the total energy

The force with respect to the position of the nuclei can be calculated similarly to the molecular case . The derivatives of all the integrals are necessary, and the derivative of the density matrix is expressed with the help of the energy-weighted density matrix. The full force is obtained as:

$`\stackrel{}{F}_{A_i}={\displaystyle \frac{E^{\mathrm{total}}}{\stackrel{}{A}_i}}=`$ (81)
$`{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \frac{T_{\mu \stackrel{}{0}\nu \stackrel{}{g}}}{\stackrel{}{A}_i}}{\displaystyle \frac{E^{\mathrm{NN}}}{\stackrel{}{A}_i}}`$ (82)
$`+{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}{\displaystyle \underset{a}{}}Z_a{\displaystyle \frac{}{\stackrel{}{A}_i}}\left[{\displaystyle \varphi _\mu (\alpha _2,\stackrel{}{r}\stackrel{}{A}_\mu ,n_2,l_2,m_2)\varphi _\nu (\alpha _1,\stackrel{}{r}\stackrel{}{A}_\nu \stackrel{}{g},n_1,l_1,m_1)A(\stackrel{}{r}\stackrel{}{A}_a)\mathrm{d}^3\mathrm{r}}\right]`$ (83)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}\{S_{\mu \stackrel{}{0}\nu \stackrel{}{g}}{\displaystyle \frac{2\pi }{3V}}{\displaystyle \underset{c}{}}{\displaystyle \underset{\stackrel{}{h},\sigma ,\tau c}{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}`$ (84)
$`{\displaystyle \frac{}{\stackrel{}{A}_i}}{\displaystyle }[\varphi _\tau (\alpha _2,\stackrel{}{r}\stackrel{}{A}_\tau ,n_2,l_2,m_2)\varphi _\sigma (\alpha _1,\stackrel{}{r}\stackrel{}{A}_\sigma \stackrel{}{h},n_1,l_1,m_1)`$ (85)
$`+{\displaystyle \underset{l=0}{\overset{L}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}{\displaystyle }\varphi _\tau (\alpha _2,\stackrel{}{r}^{}\stackrel{}{A}{}_{\tau }{}^{},n_2,l_2,m_2)\varphi _\sigma (\alpha _1,\stackrel{}{r}^{}\stackrel{}{A}_\sigma \stackrel{}{h},n_1,l_1,m_1)X_l^m(\stackrel{}{r}^{}\stackrel{}{A}_c)\mathrm{d}^3\mathrm{r}^{}\delta _l^m(\stackrel{}{A}_c,\stackrel{}{r})]r^2\mathrm{d}^3\mathrm{r}`$ (86)
$`+{\displaystyle \underset{\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}{\displaystyle \frac{C_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}}{\stackrel{}{A}_i}}`$ (87)
$`+{\displaystyle \underset{c}{}}{\displaystyle \underset{l=0}{\overset{L}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}{\displaystyle \underset{\stackrel{}{h},\tau c,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}`$ (88)
$`{\displaystyle \frac{}{\stackrel{}{A}_i}}\left[{\displaystyle }\varphi _\tau (\alpha _2,\stackrel{}{r}\stackrel{}{A}_\tau ,n_2,l_2,m_2)\varphi _\sigma (\alpha _1,\stackrel{}{r}\stackrel{}{A}_\sigma \stackrel{}{h},n_1,l_1,m_1)X_l^m(\stackrel{}{r}\stackrel{}{A_c})\mathrm{d}^3\mathrm{r}M_{l\mu \stackrel{}{0}\nu \stackrel{}{g}c}^m\right]\}`$ (89)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}^{}{\displaystyle \underset{\stackrel{}{h},\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}^{}{\displaystyle \frac{X_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}}{\stackrel{}{A}_i}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}P_{\nu \stackrel{}{g}\mu \stackrel{}{0}}^{}{\displaystyle \underset{\stackrel{}{h},\tau ,\sigma }{}}P_{\sigma \stackrel{}{h}\tau \stackrel{}{0}}^{}{\displaystyle \frac{X_{\mu \stackrel{}{0}\nu \stackrel{}{g}\tau \stackrel{}{0}\sigma \stackrel{}{h}}}{\stackrel{}{A}_i}}`$ (90)
$`+{\displaystyle \underset{\stackrel{}{g},\mu ,\nu }{}}{\displaystyle \frac{S_{\mu \stackrel{}{0}\nu \stackrel{}{g}}}{\stackrel{}{A}_i}}{\displaystyle _{BZ}}\mathrm{exp}(\mathrm{i}\stackrel{}{k}\stackrel{}{g}){\displaystyle \underset{j}{}}\{a_{\nu j}^{}(\stackrel{}{k})a_{\mu j}^{}(\stackrel{}{k})(ϵ_j^{}(\stackrel{}{k})+Q)\mathrm{\Theta }(ϵ_Fϵ_j^{}(\stackrel{}{k})Q)`$ (91)
$`+a_{\nu j}^{}(\stackrel{}{k})a_{\mu j}^{}(\stackrel{}{k})(ϵ_j^{}(\stackrel{}{k})+Q)\mathrm{\Theta }(ϵ_Fϵ_j^{}(\stackrel{}{k})Q)\}\mathrm{d}^3\mathrm{k}`$ (92)
The last addend is the energy weighted density matrix; the integral is over the first Brillouin zone.

## VI Structure of the gradient code

The present structure of the gradient code is indicated in figure 1. The first step is to compute the gradient of the Ewald energy of the nuclei in subroutine GRAMAD (the Ewald energy is computed in ENEMAD). The control module TOTGRA then first calls routines to compute the gradient of the bielectronic integrals (labeled with SHELLX$``$ as these routines will change their structure). The subroutine SHELLX$``$ calls subroutines which explicitly compute the derivatives of Coulomb and exchange integrals, and multiplies the gradients a first time with the density matrix. Back in TOTGRA again, the second multiplication with the density matrix is performed. The next step is to compute the derivatives of the multipoles (MONIRG) and to compute the energy weighted density matrix (PDIGEW). Then, the gradients of the one-electron integrals are computed (CJAT0G, CJAT1G, CJAT2G or CJAT3G, depending on the dimension). The field integrals and their gradients are now multiplied with the multipolar integrals and their gradients, and a multiplication with the density matrix is performed. This concludes the calculation of the gradients.

The structure has been simplified to focus on the most important parts. In addition, as already mentioned, the code will undergo changes during the optimization process so that a too detailed description seems to be unadvised.

## VII Examples

In tables I, II and III, we give examples of the accuracy of the gradients. First, in table I, a chain of NiO molecules is considered, with ferromagnetic ordering (all the Ni spins up) and with antiferromagnetic ordering (nearest Ni spins are antiparallel). The oxygen atoms are moved by 0.01 Å from their equilibrium positions which results in a non-vanishing force. The agreement between numerical and analytical gradient is better than 0.0001 $`\frac{E_h}{a_0}`$. As we discussed in our first article , the agreement can be improved by using stricter ”ITOL”-parameters (these are parameters which control the accuracy of the evaluation of the integrals). Indeed, when increasing these parameters, the agreement further improves up to an error of less than $`10^5\frac{E_h}{a_0}`$.

In table II, a LiF layer with a lattice constant of 5 Å is considered with one atom being displaced from its equilibrium position. The forces agree to $`2\times 10^5\frac{E_h}{a_0}`$ when default ITOL parameters (6, 6, 6, 6, 12) are used.

Finally, in table III, a three-dimensional, ferromagnetically polarized NiO solid is considered. When displacing the oxygen ions, the forces agree to better than $`2\times 10^5\frac{E_h}{a_0}`$.

As a whole, the accuracy is certainly very high and can further be improved by applying stricter cutoff (ITOL) parameters.

## VIII Conclusion

In this article, we described the implementation of analytical gradients in the code CRYSTAL. In its present form, the code is capable of computing highly accurate Hartree-Fock gradients for systems with 0,1,2 and 3-dimensional periodicity. Both closed-shell restricted Hartree-Fock as well as unrestricted Hartree-Fock calculations can be performed.

A first step of improving the efficiency of the code has been completed with the coding of gradients for the bipolar expansion, and a further enhancement of the efficiency will be one of the future directions. Of highest importance is the implementation of symmetry which will lead to high saving factors . Other targets are the implementation of gradients with respect to the lattice vector, an extension to metallic systems, and the implementation of density functional gradients.

## IX Acknowledgments

The author would like to thank CCP3 and Prof. N. M. Harrison for their interest and support of this work (EPSRC grant GR/K90661), Mr. V. R. Saunders for many helpful discussions, and Prof. R. Dovesi and the Turin group for helpful discussions and hospitality.


text:
# Optical Observations of Core-Collapse Supernovae

## INTRODUCTION

Supernovae (SNe) occur in several spectroscopically distinct varieties; see reference avf97 , for example. Type I SNe are defined by the absence of obvious hydrogen in their optical spectra, except for possible contamination from superposed H II regions. SNe II all prominently exhibit hydrogen in their spectra, yet the strength and profile of the H$`\alpha `$ line vary widely among these objects.

The early-time ($`t1`$ week past maximum brightness) spectra of SNe are illustrated in Figure 1. \[Unless otherwise noted, the optical spectra illustrated here were obtained by my group, primarily with the 3-m Shane reflector at Lick Observatory. When referring to phase of evolution, the variables $`t`$ and $`\tau `$ denote time since maximum brightness (usually in the $`B`$ passband) and time since explosion, respectively.\] The lines are broad due to the high velocities of the ejecta, and most of them have P-Cygni profiles formed by resonant scattering above the photosphere. SNe Ia are characterized by a deep absorption trough around 6150 Å produced by blueshifted Si II $`\lambda `$6355. Members of the Ib and Ic subclasses do not show this line. The presence of moderately strong optical He I lines, especially He I $`\lambda `$5876, distinguishes SNe Ib from SNe Ic.

Figure 1: Early-time spectra of SNe, showing the main subtypes.

The late-time ($`t4`$ months) optical spectra of SNe provide additional constraints on the classification scheme (Figure 2). SNe Ia show blends of dozens of Fe emission lines, mixed with some Co lines. SNe Ib and Ic, on the other hand, have relatively unblended emission lines of intermediate-mass elements such as O and Ca. At this phase, SNe II are dominated by the strong H$`\alpha `$ emission line; in other respects, most of them spectroscopically resemble SNe Ib and Ic, but with narrower emission lines. The late-time spectra of SNe II show substantial heterogeneity, as do the early-time spectra.

To a first approximation, the light curves of SNe I are all broadly similar lei91a . SNe Ib usually have slower decline rates than SNe Ic; however, SNe Ic may come in “slow” and “fast” varieties clo97a ; clo97b . The light curves of SNe II exhibit much dispersion pat93 , though it is useful to subdivide the majority of them into two relatively distinct subclasses bar79 ; dog85 . The light curves of SNe II-L (“linear”) generally resemble those of SNe I, with a steep decline after maximum brightness followed by a slower exponential tail. In contrast, SNe II-P (“plateau”) remain within $`1`$ mag of maximum brightness for an extended period. The light curve of SN 1987A, albeit atypical, was generically related to those of SNe II-P.

Figure 2: Late-time spectra of SNe. At even later phases, SN 1987A was dominated by strong emission lines of H$`\alpha `$, \[O I\], \[Ca II\], and the Ca II near-infrared triplet.

The locations at which SNe occur provide important clues to their nature, and to the mass of their progenitor stars. SNe II, Ib, and Ic have never been seen in elliptical galaxies, and rarely if ever in S0 galaxies. They are generally in or near spiral arms and H II regions van96 , implying that their progenitors must have started their lives as massive stars ($`10M_{}`$). The progenitors of SNe II are thought to suffer core collapse and subsequently “rebound” with help from neutrinos arn89 ; bur00 , leaving a neutron star or perhaps in some cases a black hole bro94 . Most workers now believe that SNe Ib/Ic are produced by the same mechanism as SNe II, except that the progenitors were stripped of their hydrogen (SN Ib) and possibly helium (SN Ic) envelopes prior to exploding, either via mass transfer to companion stars nom94 ; woo95  or through winds (e.g., woo93 ; swa93 ). White dwarf models have been discussed bra86  but are unlikely.

Figure 3: Montage of spectra of SN 1992H in NGC 5377. Epochs (days) are given relative to the estimated time of explosion, February 8, 1992.

## SUBCLASSES OF TYPE II SUPERNOVAE

Most SNe II-P seem to have a relatively well-defined spectral development, as shown in Figure 3 for SN 1992H (see also reference clo96 ). At early times the spectrum is nearly featureless and very blue, indicating a high color temperature ($``$ 10,000 K). He I $`\lambda `$5876 with a P-Cygni profile is sometimes visible. The temperature rapidly decreases with time, reaching $`5000`$ K after a few weeks, as expected from the adiabatic expansion and associated cooling of the ejecta. It remains roughly constant at this value during the plateau (the photospheric phase), while the hydrogen recombination wave moves through the massive ($`10M_{}`$) hydrogen ejecta and releases the energy deposited by the shock. At this stage strong Balmer lines and Ca II H&K with well-developed P-Cygni profiles appear, as do weaker lines of Fe II, Sc II, and other iron-group elements. The spectrum gradually takes on a nebular appearance as the light curve drops to the late-time tail; the continuum fades, but H$`\alpha `$ becomes very strong, and prominent emission lines of \[O I\], \[Ca II\], and Ca II also appear.

Few SNe II-L have been observed in as much detail as SNe II-P. Figure 4 shows the spectral development of SN 1979C bra81 , an unusually luminous member of this subclass. Near maximum brightness the spectrum is very blue and almost featureless, with a slight hint of H$`\alpha `$ emission. A week later, H$`\alpha `$ emission is more easily discernible, and low-contrast P-Cygni profiles of Na I, H$`\beta `$, and Fe II have appeared. By $`t1`$ month, the H$`\alpha `$ emission line is very strong but still devoid of an absorption component, while the other features clearly have P-Cygni profiles. Strong, broad H$`\alpha `$ emission dominates the spectrum at $`t7`$ months, and \[O I\] $`\lambda \lambda `$6300, 6364 emission is also present. Several authors whe90 ; avf91a ; sch96  have speculated that the absence of H$`\alpha `$ absorption spectroscopically differentiates SNe II-L from SNe II-P, but the small size of the sample of well-observed objects precluded definitive conclusions.

Figure 4: Montage of spectra of SN 1979C in NGC 4321, from reference bra81 ; reproduced with permission. Epochs (days) are given relative to the date of maximum brightness, April 15, 1979.

The progenitors of SNe II-L are generally believed to have relatively low-mass hydrogen envelopes (a few $`M_{}`$); otherwise, they would exhibit distinct plateaus, as do SNe II-P. On the other hand, they may have more circumstellar gas than do SNe II-P, and this could give rise to the emission-line dominated spectra. They are often radio sources sra90 ; moreover, the ultraviolet excess (at $`\lambda 1600`$ Å) seen in SNe 1979C and 1980K may be produced by inverse Compton scattering of photospheric radiation by high-speed electrons in shock-heated ($`T10^9`$ K) circumstellar material fra82 ; fra84 . Finally, the light curves of some SNe II-L reveal an extra source of energy: after declining exponentially for several years, the H$`\alpha `$ flux of SN 1980K reached a steady level, showing little if any decline thereafter uom86 ; lei91b . The excess almost certainly comes from the kinetic energy of the ejecta being thermalized and radiated due to an interaction with circumstellar matter che90 ; lei94 .

The very late-time optical recovery of SNe 1979C and 1980K lei91b ; fes95 ; fes99  and other SNe II-L supports the idea of ejecta interacting with circumstellar material. The spectra consist of a few strong, broad emission lines such as H$`\alpha `$, \[O I\] $`\lambda \lambda `$6300, 6364, and \[O III\] $`\lambda \lambda `$4959, 5007. A Hubble Space Telescope (HST) ultraviolet spectrum of SN 1979C reveals some prominent, double-peaked emission lines with the blue peak substantially stronger than the red, suggesting dust extinction within the expanding ejecta fes99 . The data show general agreement with the emission lines expected from circumstellar interaction che94 , but the specific models that are available show several differences with the observations. For example, we find higher electron densities ($`10^5`$ to $`10^7`$ cm-3), resulting in stronger collisional de-excitation than assumed in the models. These differences can be used to further constrain the nature of the progenitor star. Note that based on photometry of the stellar populations in the environment of SN 1979C (from HST images), the progenitor of the SN was at most 10 million years years old, so its initial mass was probably 17–18 $`M_{}`$ van99 .

During the past decade, there has been the gradual emergence of a new, distinct subclass of SNe II avf91a ; avf91b ; sch90 ; lei94  whose ejecta are believed to be strongly interacting with dense circumstellar gas, even at early times (unlike SNe II-L). The derived mass-loss rates for the progenitors can exceed $`10^4M_{}`$ yr-1 chu94 . In these objects, the broad absorption components of all lines are weak or absent throughout their evolution. Instead, their spectra are dominated by strong emission lines, most notably H$`\alpha `$, having a complex but relatively narrow profile. Although the details differ among objects, H$`\alpha `$ typically exhibits a very narrow component (FWHM $`200`$ km s-1) superposed on a base of intermediate width (FWHM $``$ 1000–2000 km s-1; sometimes a very broad component (FWHM $``$ 5000–10,000 km s-1) is also present. This subclass was christened “Type IIn” sch90 , the “n” denoting “narrow” to emphasize the presence of the intermediate-width or very narrow emission components. Representative spectra of five SNe IIn are shown in Figure 5, with two epochs for SN 1994Y.

The early-time continua of SNe IIn tend to be bluer than normal. Occasionally He I emission lines are present in the first few spectra (e.g., SN 1994Y in Figure 5). Very narrow Balmer absorption lines are visible in the early-time spectra of some of these objects, often with corresponding Fe II, Ca II, O I, or Na I absorption as well (e.g., SNe 1994W and 1994ak in Figure 5). Some of them are unusually luminous at maximum brightness, and they generally fade quite slowly, at least at early times. The equivalent width of the intermediate H$`\alpha `$ component can grow to astoundingly high values at late times. The great diversity in the observed characteristics of SNe IIn provides clues to the various degrees and forms of mass loss late in the lives of massive stars.

Figure 5: Montage of spectra of SNe IIn. Epochs are given relative to the estimated dates of explosion.

## TYPE II SUPERNOVA IMPOSTORS?

The peculiar SN IIn 1961V (“Type V” according to Zwicky zwi65 ) had probably the most bizarre light curve ever recorded. (SN 1954J, also known as “Variable 12” in NGC 2403, was similar hum94 .) Its progenitor was a very luminous star, visible in many photographs of the host galaxy (NGC 1058) prior to the explosion. Perhaps SN 1961V was not a genuine supernova (defined to be the violent destruction of a star at the end of its life), but rather the super-outburst of a luminous blue variable such as $`\eta `$ Carinae goo89 ; avf95 .

A related object may have been SN IIn 1997bs, the first SN discovered in the Lick Observatory Supernova Search (LOSS, described later in this review). Its spectrum was peculiar (Figure 6), consisting of narrow Balmer and Fe II emission lines superposed on a featureless continuum. Its progenitor was discovered in an HST archival image of the host galaxy van00 . It is a very luminous star ($`M_V7.4`$ mag), and it didn’t brighten as much as expected for a SN explosion ($`M_V13`$ at maximum). These data suggest that SN 1997bs may have been like SN 1961V — that is, a supernova impostor. The real test will be whether the star is still visible in future HST images obtained years after the outburst.

Figure 6: Spectrum of SN 1997bs, obtained on April 16, 1997 UT.

## TYPE Ib AND Ic SUPERNOVAE

In the 1960s, Bertola and collaborators ber64 ; ber65  recognized that not all SNe I are of the “classical” variety (now known as SNe Ia) with a strong absorption trough near 6150 Å. By the mid-1980s these objects came to be known as SNe Ib eli85 , and He I lines were identified in their early-time spectra har87 . Gradually it became clear that SNe Ib constitute a heterogeneous subclass, with substantial variations in the observed He I strengths in spectra obtained around maximum brightness. Wheeler & Harkness (whe86 ; see also har87 ) suggested that SNe Ib should actually be divided into two separate categories: SNe Ib are those showing strong He I absorption lines (especially He I $`\lambda `$5876) in their early-time photospheric spectra, while SNe Ic are those in which He I is not easily discernible. However, they modeled SNe Ic in the same physical way as SNe Ib whe87 , but with different relative concentrations of He and O in the envelope.

A large, comprehensive study of SNe Ib and SNe Ic was recently completed by my group mat01 . The relative depths of the helium absorption lines in the spectra of the SNe Ib appear to provide a measurement of the temporal evolution of the supernova, with He I $`\lambda `$5876 and He I $`\lambda `$7065 growing in strength relative to He I $`\lambda `$6678 over time. Some SNe Ic show evidence for weak He I absorption, but most do not. Aside from the presence or absence of the helium lines, there are other spectroscopic differences between SNe Ib and SNe Ic. On average, the O I $`\lambda `$7774 line is stronger in SNe Ic than in SNe Ib. In addition, the SNe Ic have distinctly broader emission lines at late times, indicating either a consistently larger explosion energy and/or a lower envelope mass for SNe Ic than for SNe Ib. These results are consistent with the idea that the progenitors of SNe Ic are massive stars that have lost more of their envelope (i.e., much of the helium layer) than the progenitors of SNe Ib. The general hypothesis that SNe Ib/Ic have “stripped” progenitors is greatly supported by the discovery of links between SNe II and SNe Ib/Ic, as I discuss next.

## LINKS BETWEEN TYPE II AND TYPE Ib/Ic SUPERNOVAE

Filippenko avf88  presented the case of SN 1987K, which appeared to be a link between SNe II and SNe Ib. Near maximum brightness, it was undoubtedly a SN II, but with rather weak photospheric Balmer and Ca II lines. Many months after maximum brightness, its spectrum was essentially that of a SN Ib. The simplest interpretation is that SN 1987K had a meager hydrogen atmosphere at the time it exploded; it would naturally masquerade as a SN II for a while, and as the expanding ejecta thinned out the spectrum would become dominated by emission from deeper and denser layers. The progenitor was probably a star that, prior to exploding via iron core collapse, lost almost all of its hydrogen envelope either through mass transfer onto a companion or as a result of stellar winds. Such SNe were dubbed “SNe IIb” by Woosley et al. woo87 , who had proposed a similar preliminary model for SN 1987A before it was known to have a massive hydrogen envelope.

The data for SN 1987K (especially its light curve) were rather sparse, making it difficult to model in detail. Fortunately, the Type II SN 1993Jin NGC 3031 (M81) came to the rescue, and was studied in greater detail than any supernova since SN 1987A whe96 . Its light curves ric96  and spectra avf93 ; avf94 ; mat00a ; mat00b  amply supported the hypothesis that the progenitor of SN 1993J probably had a low-mass (0.1–0.6 $`M_{}`$) hydrogen envelope above a $`4M_{}`$ He core nom93 ; pod93 ; woo94 . Figure 7 shows several early-time spectra of SN 1993J, showing the emergence of He I features typical of SNe Ib. Considerably later (Figure 8), the H$`\alpha `$ emission nearly disappeared, and the spectral resemblance to SNe Ib was strong. The general consensus is that its initial mass was $`15M_{}`$. A star of such low mass cannot shed nearly its entire hydrogen envelope without the assistance of a companion star. Thus, the progenitor of SN 1993J probably lost most of its hydrogen through mass transfer to a bound companion 3–20 AU away. In addition, part of the gas may have been lost from the system. Had the progenitor lost essentially all of its hydrogen prior to exploding, it would have had the optical characteristics of SNe Ib. There is now little doubt that most SNe Ib, and probably SNe Ic as well, result from core collapse in stripped, massive stars, rather than from the thermonuclear runaway of white dwarfs.

Figure 7: Early-time spectral evolution of SN 1993J. A comparison with the Type Ib SN 1984L is shown at bottom, demonstrating the presence of He I lines in SN 1993J. The explosion date was March 27.5, 1993.

SN 1993J held several more surprises. Observations at radio van94  and X-ray suz95  wavelengths revealed that the ejecta are interacting with relatively dense circumstellar material fra96 , probably ejected from the system during the course of its pre-SN evolution. Optical evidence for this interaction also began emerging at $`\tau 10`$ months: the H$`\alpha `$ emission line grew in relative prominence, and by $`\tau 14`$ months it had become the dominant line in the spectrum avf94 ; pat95 ; fin95 , consistent with models che94 . Its profile was very broad (FWHM $``$ 17,000 km s-1; Figure 8) and had a relatively flat top, but with prominent peaks and valleys whose likely origin is Rayleigh-Taylor instabilities in the cool, dense shell of gas behind the reverse shock che92 . Radio VLBI measurements show that the ejecta are circularly symmetric, but with significant emission asymmetries mar95 , possibly consistent with the asymmetric H$`\alpha `$ profile seen in some of the spectra avf94 .

Figure 8: In the top spectrum, which shows SN 1993J about 7 months after the explosion, H$`\alpha `$ emission is very weak; the resemblance to spectra of SNe Ib is striking. A year later (bottom), however, H$`\alpha `$ was once again the dominant feature in the spectrum (which was scaled for display purposes).

## SPECTROPOLARIMETRY OF SUPERNOVAE

Spectropolarimetry of SNe can be used to probe their geometry. The basic question is whether SNe are round, and the idea is simple: A hot young SN atmosphere is dominated by electron scattering, which by its nature is highly polarizing. For an unresolved, spherical source, however, the directional components of the electric vectors cancel exactly, yielding zero net linear polarization. Conversely, if the source is aspherical, incomplete cancellation occurs and a net linear polarization results, with typical continuum polarizations of $`1\%`$ expected for objects with moderate ($`20\%`$) asphericity (e.g., hof91 ), although the exact amount of polarization observed is also sensitive to the viewing angle. Spectropolarimetry is important for a full understanding of the physics of SN explosions and can also provide information on the circumstellar environment of SNe.

My group obtained spectropolarimetry of one object from each of the major SN types and subtypes leo00a , generally with the Keck-II 10-m telescope (but in a few important cases with the Lick 3-m Shane reflector). Most of the objects exhibit a change in both the magnitude and direction of the polarization across strong lines, especially the absorption troughs of the strongest P-Cygni lines. This may result from global asymmetry of the electron-scattering atmosphere and/or the underlying continuum region mcc84 ; hof91 .

We have studied SN IIn 1998S in some detail leo00b  (see also wan01 ); its optical spectrum is dominated by strong, multi-component emission lines, thought to be produced by an intense interaction between the supernova and its dense circumstellar environment. We combined our early-time (3 days after discovery) spectropolarimetric observation with total flux spectra spanning nearly 500 days. We measure an intrinsic continuum polarization of $`p3\%`$ (one of the highest yet found for a SN), suggesting a global asphericity of $`45\%`$ from the models of Höflich hof91 . The line profiles favor a ring-like geometry for the circumstellar gas, generically similar to what is seen directly in SN 1987A, although much denser and closer to the progenitor in SN 1998S.

We also found that one month after exploding, the Type IIb SN 1993J had a polarized flux spectrum resembling spectra of SNe Ib, with prominent He I lines tra97 . The data are consistent with models in which the polarization is produced by an asymmetric He core configuration of material. It is interesting that the percent polarization may increase with decreasing envelope mass, along the sequence Type II, IIb, Ib, and Ic whe00 ; wan01 ; leo01b , suggesting that asymmetries in massive stars become more pronounced as one probes deeper into the core.

SN 1999em, an extremely bright ($`m_V13.5`$) SN II-P, provided the rare opportunity to study the geometry of a “normal” core-collapse event at multiple epochs leo01a . We obtained spectropolarimetry at 3 plateau-phase epochs, and then one final epoch long after it had dropped off the plateau (Fig. 9). A very low but temporally increasing polarization level suggests a substantially spherical geometry at early times that becomes more aspherical at late times as ever-deeper layers of the ejecta are revealed. We speculate that the thick hydrogen envelope intact at the time of explosion in SNe II-P might serve to dampen the effects of an intrinsically aspherical explosion. The increase in asphericity seen at later times is consistent with the trend identified above among stripped-envelope core-collapse SNe: the deeper we peer, the more evidence we find for asphericity. The natural conclusion that it is explosion asymmetry that is responsible for the polarization has fueled the idea that some core-collapse SNe produce gamma-ray bursts (GRB; e.g., blo99 ) through the action of a jet of material aimed fortuitously at the observer, the result of a “bipolar,” jet-induced, SN explosion kho99 ; whe00 .

Figure 9: Montage of the observed optical polarization of SN 1999em, with arbitrarily scaled total flux spectra (thick lines) overplotted for comparison of features. The polarization spectra are binned 5 Å bin-1 to improve the S/N ratio. See reference leo01b  for details.

## SUPERNOVAE ASSOCIATED WITH GAMMA-RAY BURSTS?

At least a small fraction of gamma-ray bursts (GRBs) may be associated with nearby SNe. Probably the most compelling example thus far is that of SN 1998bw and GRB 980425 (e.g., gal98 ; iwa98 ; woo99 ; sta00 ), which were temporally and spatially coincident. SN 1998bw was, in many ways, an extraordinary SN; it was very luminous at optical and radio wavelengths, and it showed evidence for relativistic outflow. Its bizarre optical spectrum is often classified as that of a SN Ic, but the object should be called a “peculiar SN Ic” if not a subclass of its own; the spectrum was distinctly different from that of a normal SN Ic.

Models suggest that SNe associated with GRBs are highly asymmetric; thus, spectropolarimetry should provide some useful tests. In particular, perhaps objects such as SN 1998S, mentioned above, would have been seen as GRBs had their rotation axis been pointed in our direction. That of SN 1998S was almost certainly not aligned with us leo00b ; both the spectropolarimetry and the appearance of double-peaked H$`\alpha `$ emission suggest an inclined view, rather than pole-on.

Figure 10: Spectral evolution of SN 1999E, which may have been associated with GRB 980910.

The case of GRB 970514 and the very luminous SN IIn 1997cy is also interesting ger00 ; tur00 ; there is a reasonable possibility that the two objects were associated. The optical spectrum of SN 1997cy was highly unusual, and bore some resemblance to that of SN 1998bw, though there were some differences as well. SN 1999E, which might be linked with GRB 980910 but with large uncertainties tho99 , also had an optical spectrum similar to that of SN 1997cy avf99 ; cap99a ; see Figure 10. The undulations are very broad, indicating high ejection velocities. Besides H$`\alpha `$, secure line identifications are difficult, though some of the emission features seem to be associated with oxygen and calcium. Perhaps SN 1999E was produced by the highly asymmetric collapse of a carbon-oxygen core.

## SUPERNOVA RATES

I conclude with a brief summary of supernova rates. The derived rates of various types of SNe as a function of Hubble type are still quite uncertain; no single search has found a sufficiently large sample of SNe, and combining different searches can be misleading due to differences in the dominant selection effects. Probably the best existing determination of supernova rates is that of Capellaro et al. cap99b , who combined the results of five searches while being attentive to various selection effects. Four of the searches had been done with photographic plates, while the fifth was a visual search by the Rev. Robert Evans. Given that photographic searches can be rather insensitive to SNe near the central regions of galaxies, and that the Evans visual search did not go very deep, we must avoid the temptation to take the overall results (see table below) too literally. Nevertheless, they can be used as a rough guideline.

SNu \[#(100 yr)$`{}_{}{}^{1}(10^{10}L_{}^B)_{}^{1}`$\]

| Galaxy type | Ia | Ib/c | II | All |
| --- | --- | --- | --- | --- |
| E-S0 | $`0.18\pm 0.06`$ | $`<0.01`$ | $`<0.02`$ | $`0.18\pm 0.06`$ |
| S0a-Sb | $`0.18\pm 0.07`$ | $`0.11\pm 0.06`$ | $`0.42\pm 0.19`$ | $`0.72\pm 0.21`$ |
| Sbc-Sd | $`0.21\pm 0.08`$ | $`0.14\pm 0.07`$ | $`0.86\pm 0.35`$ | $`1.21\pm 0.37`$ |
| Others | $`0.40\pm 0.16`$ | $`0.22\pm 0.16`$ | $`0.65\pm 0.39`$ | $`1.26\pm 0.45`$ |
| All | $`0.20\pm 0.06`$ | $`0.08\pm 0.04`$ | $`0.40\pm 0.19`$ | $`0.68\pm 0.20`$ |

*Note:* $`H_0=75`$ km s-1 Mpc-1; scale by ($`H_0/75)^2`$.

Others include Sm, Irr, Pec.

My group is trying to remedy the situation by conducting a long-term search for nearby SNe (with redshifts generally less than 5000 km s-1) in a uniform manner wli00 ; avf01 . Special emphasis is placed on finding them well before maximum brightness. We are using the Katzman Automatic Imaging Telescope (KAIT) at Lick Observatory, a fully robotic 0.75-m reflector equipped with a CCD imaging camera. Its telescope control system checks the weather, opens the dome, points to the desired objects, finds and acquires guide stars (only for long integrations), exposes, stores the data, and manipulates the data without human intervention. There is a filter wheel with 20 slots, including UBVRI. Five-minute guided exposures yield $`R20`$ mag.

A limit of about 19 mag ($`4\sigma `$) is reached in the 25-second unfiltered, unguided exposures of our Lick Observatory Supernova Search (LOSS). We observe up to 1200 galaxies per night, and try to cycle back to the same galaxies after 3–4 nights. Our software automatically subtracts new images from old ones and identifies supernova candidates that are subsequently examined by undergraduate research assistants. LOSS found 20 SNe in 1998, 40 in 1999, and 36 in 2000, making KAIT the world’s most successful search engine for nearby SNe. (Note that we found SN 2000A and SN 2001A — and hence the first supernova of the new millennium regardless of one’s definition of the turn of the millennium!) After a few more years, we hope to have found enough SNe to provide a meaningful update on the supernova rates given above.

Multi-filter follow-up photometry is conducted of the most important SNe, and all objects are monitored in unfiltered mode. A Web page describing LOSS is at http://astro.berkeley.edu/$``$bait/kait.html . As a byproduct of LOSS, we also find novae in the Local Group, comets, asteroids, and cataclysmic variables.

## ACKNOWLEDGMENTS

My recent research on SNe has been financed by the NSF, most recently through grant AST-9987438, as well as by NASA grants GO-7821, GO-8243, GO-8648, AR-6371, and AR-8006 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA Contract NAS5-26555. Additional funding was provided by NASA/Chandra grant GO-0-1009C. KAIT and its associated science have been made possible with funding or donations from NSF, NASA, the Sylvia and Jim Katzman Foundation, Sun Microsystems Inc., Lick Observatory, the Hewlett-Packard Company, Photometrics Ltd., AutoScope Corporation, and the University of California. This paper was written while I was supported by a Guggenheim Fellowship. I am grateful to the students and postdocs who have worked with me on SNe over the past 15 years for their assistance and discussions.