text:
# Measurement-induced spin squeezing in a cavity

## I Introduction

Squeezed spin states (SSS) Kitagawa , in which quantum uncertainty of the collective spin in a direction orthogonal to the mean spin vector is suppressed below the standard quantum limit, have recently been realized experimentally Hald ; Kuzmich00 ; Juls . In the SSS of atomic systems, entanglement between atoms is crucial for suppressing fluctuations in collective atomic spin, and therefore such systems are of potential interest in quantum information processing as well as in atomic interferometry.

A number of methods to generate SSS in atomic systems have been proposed, and these can be classified into several categories. A few examples that are relevant to the subject of the present paper are: (i) Squeezing is transferred between photons and atoms via their interactions SaitoL ; Kuzmich97 ; SaitoA ; Poulsen . Using this property, SSS can be realized in an atomic system by copying the squeezing from the system of photons to the atoms. This scheme has been demonstrated experimentally Hald . (ii) A probe laser beam propagating through atoms establishes entanglement between the atoms, with measurement of the probe laser inducing state reduction in the atomic states, resulting in SSS Kuzmich98 ; Takahashi ; Duan . This kind of measurement-induced-squeezing scheme has also been employed in experiments Kuzmich00 ; Juls . (iii) The interaction between atoms and photons confined in a high-Q cavity entangles them, and states of both atoms and photons develop into squeezed states Wineland ; Waka ; SaitoL ; SaitoA .

In the method (iii) proposed by Wineland et al. Wineland  and by us Waka ; SaitoL ; SaitoA , photons prepared in a coherent state and atoms prepared in a coherent spin state undergo interactions in a high-Q cavity, and after a certain interaction time SSS is realized temporarily, but not in a steady state. Precise control over the initial state and interaction time is therefore required. In the present paper, a novel scheme to generate SSS is proposed that combines methods (iii) and (ii) to circumvent the stringent requirements of the precise control needed for method (iii). SSS is thus obtained as a “quasi-steady state”, regardless of the initial state, by performing continuous measurement of photons in the cavity. This drastically reduces experimental difficulties in producing SSS. Moreover, a degree of spin squeezing proportional to $`N^{1/3}`$, where $`N`$ is the number of atoms, can be achieved using this method, in contrast to the $`N^{1/4}`$ attainable by method (iii) SaitoA .

This paper is organized as follows. Section II describes the system we consider and develops a theory of continuous measurement of photon number in the cavity that is being pumped and contains atoms to be squeezed. Section III shows that the atomic state in the cavity develops into an SSS by continuous measurement of photons in the cavity. Section IV discusses possible experimental schemes and concludes this paper.

## II Continuous Measurement of Photon Number in a Pumped Cavity with Atoms

### II.1 Description of the system

An operational theory of continuous photon counting has been proposed by Srinivas and Davies Srinivas , and microscopic modeling of this theory has been reported in Ref. Imoto . The nonunitary state evolution due to measurement back action associated with the Srinivas-Davies model was investigated in Ref. Ueda . In these studies, continuous measurement of single-mode photons in a high-Q cavity was considered. We extend this model to the case of a pumped cavity containing $`N`$ two-level atoms as illustrated in Fig. 1.

We first consider the system without photon counting, namely, without the probe atoms and detector shown in Fig. 1. The high-Q cavity is assumed to sustain single-mode photons of frequency $`\omega _\mathrm{c}`$ and be pumped by a coherent light source with an effective amplitude $`A`$ and frequency $`\omega _\mathrm{p}`$ inside the cavity. The cavity contains two-level atoms with transition frequency $`\omega _\mathrm{a}`$, which couple to photons via the Jaynes-Cummings interaction JCM  with coupling constant $`g`$. We assume that $`g`$ is constant and the same for all the atoms, that is, we consider the Lamb-Dicke limit, where all atoms experience the same photon field. It is then convenient to introduce a set of collective spin operators defined by

$`\widehat{S}_z`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{N}{}}}(|\mathrm{e}_j\mathrm{e}_j||\mathrm{g}_j\mathrm{g}_j|),`$ (1)
$`\widehat{S}_+`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{N}{}}}|\mathrm{e}_j\mathrm{g}_j|\widehat{S}_{}^{},`$ (2)
where $`N`$ is the number of atoms, and $`|\mathrm{g}_j`$ and $`|\mathrm{e}_j`$ are the ground and excited states of the $`j`$th atom. We shall also use other collective spin operators defined by $`\widehat{S}_x(\widehat{S}_++\widehat{S}_{})/2`$ and $`\widehat{S}_y(\widehat{S}_+\widehat{S}_{})/(2i)`$.

The Hamiltonian for the system is given by

$$\widehat{H}_0/\mathrm{}=\omega _\mathrm{c}\widehat{n}+\omega _\mathrm{a}\widehat{S}_z+A(e^{i\omega _\mathrm{p}t}\widehat{a}+e^{i\omega _\mathrm{p}t}\widehat{a}^{})+g(\widehat{a}\widehat{S}_++\widehat{a}^{}\widehat{S}_{}),$$
(3)
where $`\widehat{a}`$ ($`\widehat{a}^{}`$) is the annihilation (creation) operator of the photon field, $`\widehat{n}\widehat{a}^{}\widehat{a}`$ is the photon-number operator, and $`A`$ and $`g`$ are taken to be real and positive without loss of generality. In Eq. (3), the first and the second terms on the right-hand side (rhs) describe the free Hamiltonian for the photon field and the collective atomic system, respectively; the third term describes the effect of coherent pumping of the photon field, and the last term describes the interaction between the photon field and the atomic system. In the rotating frame of reference $`e^{i\omega _\mathrm{p}(\widehat{n}+\widehat{S}_z)}|\psi `$, the Hamiltonian (3) is time-independent:

$$\widehat{H}_0^{\mathrm{rot}}/\mathrm{}=\delta _{\mathrm{cp}}\widehat{n}+\delta _{\mathrm{ap}}\widehat{S}_z+A(\widehat{a}+\widehat{a}^{})+g(\widehat{a}\widehat{S}_++\widehat{a}^{}\widehat{S}_{}),$$
(4)
where $`\delta _{\mathrm{cp}}\omega _\mathrm{c}\omega _\mathrm{p}`$ and $`\delta _{\mathrm{ap}}\omega _\mathrm{a}\omega _\mathrm{p}`$.

We now consider continuous photon counting for the system described by the Hamiltonian (4). Experimentally, photon counting may be performed using a photodetector placed in the immediate vicinity of the cavity, where photons leaking out of the cavity are detected. This kind of continuous photodetection can be simulated using a sequence of probe atoms that are injected through the cavity Imoto , as illustrated in Fig. 1. The probe atoms are prepared in the ground state and injected into the cavity, with internal states detected upon leaving the cavity. If a probe atom is found to be in an excited state, we can say that one photon has been detected (one-count process). If it is found to have remained in the ground state, we can say that no photon has been detected (no-count processImoto ; Ueda . In the remaining part of this section we discuss nonunitary state evolution of the system in the one-count process, the no-count process, and more general processes.

### II.2 State evolution in the one-count process

If the passage time of each probe atom through the cavity is much shorter than the time scale of the state evolution described by Eq. (4), the density operator of the entire system (i.e., cavity photons plus atoms excluding the probe atoms) immediately after the one-count process at time $`t`$ is given by (see Appendix A for the derivation)

$$\widehat{\rho }(t_+)=\frac{\widehat{a}\widehat{\rho }(t)\widehat{a}^{}}{\mathrm{Tr}[\widehat{a}\widehat{\rho }(t)\widehat{a}^{}]},$$
(5)
where $`t_+`$ denotes a time infinitesimally later than $`t`$ and $`\mathrm{Tr}`$ the trace of the entire system. The denominator on the rhs of Eq. (5) multiplied by $`\gamma \mathrm{\Delta }t`$ gives the probability that a photon is detected between $`t`$ and $`t+\mathrm{\Delta }t`$, where $`\gamma `$ is the photodetection rate defined by Eq. (37). The expression (5) takes the same form as that of the Srinivas-Davies model Srinivas . Immediately after the one-count process (5), the expectation value of an arbitrary operator $`\widehat{𝒪}`$ is given by

$$\mathrm{Tr}[\rho (t_+)\widehat{𝒪}]\widehat{𝒪}_{t_+}=\frac{\widehat{a}^{}\widehat{𝒪}\widehat{a}_t}{\widehat{n}_t}.$$
(6)
When $`\widehat{𝒪}`$ is the photon-number operator $`\widehat{n}`$, we obtain

$$\widehat{n}_{t_+}=\widehat{n}_t1+\frac{(\mathrm{\Delta }\widehat{n})^2_t}{\widehat{n}_t},$$
(7)
where $`\mathrm{\Delta }\widehat{n}\widehat{n}\widehat{n}`$. The one-count process thus decreases the mean photon number by one and increases it by the Fano factor $`(\mathrm{\Delta }\widehat{n})^2_t/\widehat{n}_t`$ of the photon field immediately before the one-count process Ueda .

When $`\widehat{𝒪}`$ is some spin operator $`\widehat{𝒮}`$, which commutes with $`\widehat{a}`$, Eq. (6) becomes

$$\widehat{𝒮}_{t_+}=\frac{\widehat{𝒮}\widehat{n}_t}{\widehat{n}_t}=\widehat{𝒮}_t+\frac{\mathrm{\Delta }\widehat{𝒮}\mathrm{\Delta }\widehat{n}_t}{\widehat{n}_t}.$$
(8)
This indicates that upon the detection of a photon, the expectation value of the collective spin operator changes by an amount that depends on the correlation between the photons and atoms. As an example, let us consider an entangled state between an atom and a photon, where the ground and excited states of the atom are denoted by $`|\mathrm{g}`$ and $`|\mathrm{e}`$, and the vacuum and single-photon state by $`|0`$ and $`|1`$. Suppose that the state of the system before measurement is given by $`(|\mathrm{g}|1+|\mathrm{e}|0)/\sqrt{2}`$. Then $`\widehat{S}_z`$ decreases from $`1/2`$ to $`0`$ by the one-count process. This is because the measurement reveals that the cavity had contained a photon, and hence that the state of the system after the one-count process is given by $`|\mathrm{g}|0`$. Similarly, $`\widehat{S}_z`$ increases upon the one-count process from $`1/2`$ to $`1`$ for an entangled initial state $`(|\mathrm{g}|0+|\mathrm{e}|1)/\sqrt{2}`$. In contrast, when the initial state is separable, i.e., $`\widehat{𝒮}\widehat{n}_t=\widehat{𝒮}_t\widehat{n}_t`$, the atomic quantities do not change by the one-count process.

### II.3 State evolution in the no-count process

When no photon is detected between $`t`$ and $`t+T`$, the density operator is shown to evolve according to the following equation (see Appendix A for the derivation),

$$\widehat{\rho }(t+T)=\frac{e^{i\widehat{H}_{\mathrm{nc}}T/\mathrm{}}\widehat{\rho }(t)e^{i\widehat{H}_{\mathrm{nc}}^{}T/\mathrm{}}}{\mathrm{Tr}[e^{i\widehat{H}_{\mathrm{nc}}T/\mathrm{}}\widehat{\rho }(t)e^{i\widehat{H}_{\mathrm{nc}}^{}T/\mathrm{}}]},$$
(9)
where $`\widehat{H}_{\mathrm{nc}}\widehat{H}_0^{\mathrm{rot}}i\mathrm{}\gamma \widehat{n}/2`$ is the non-Hermitian operator that governs the state evolution during the no-count process. The denominator on the rhs of Eq. (9) gives the probability that no photon is detected between $`t`$ and $`t+T`$. It follows from Eq. (9) that the master equation for the no-count process is given by

$$\frac{\widehat{\rho }(t)}{t}=\frac{i}{\mathrm{}}[\widehat{\rho }(t)\widehat{H}_{\mathrm{nc}}^{}\widehat{H}_{\mathrm{nc}}\widehat{\rho }(t)]+\gamma \widehat{n}_t\widehat{\rho }(t).$$
(10)
In fact, it can be shown by direct substitution of Eq. (9) into Eq. (10) that Eq. (9) is the solution of Eq. (10).

### II.4 State evolution in the referring and non-referring measurement processes

Time evolution under continuous photodetection is thus a stochastic process, in which the state evolution associated with one-count and no-count processes are described by Eqs. (5) and (9), respectively. The probability of the one-count process occurring between $`t`$ and $`t+dt`$ is $`\gamma \widehat{n}_tdt`$ and that of the no-count process is $`1\gamma \widehat{n}_tdt`$. The number of photons detected and the times at which they are detected, therefore, differ run by run even for the same initial state. The measurement process in which we read the outcome of the measurement at every instant of time is referred to as referring measurement process Ueda . When the photons are detected at $`t_1,t_2,\mathrm{},t_n`$ from $`t=0`$ to $`t=T`$, the time evolution of the density operator is given by

$$\widehat{\rho }(T)=\frac{\widehat{P}(T;t_1,\mathrm{},t_n)\widehat{\rho }(0)\widehat{P}^{}(T;t_1,\mathrm{},t_n)}{\mathrm{Tr}[\widehat{P}(T;t_1,\mathrm{},t_n)\widehat{\rho }(0)\widehat{P}^{}(T;t_1,\mathrm{},t_n)]},$$
(11)
where

$`\widehat{P}(T;t_1,\mathrm{},t_n)`$ $``$ $`e^{i\widehat{H}_{\mathrm{nc}}(Tt_n)/\mathrm{}}\widehat{a}e^{i\widehat{H}_{\mathrm{nc}}(t_nt_{n1})/\mathrm{}}\widehat{a}`$ (12)
$`\mathrm{}\widehat{a}e^{i\widehat{H}_{\mathrm{nc}}(t_2t_1)/\mathrm{}}\widehat{a}e^{i\widehat{H}_{\mathrm{nc}}t_1/\mathrm{}}.`$
If the detector is switched on but we do not read out the measurement results, the state evolves into a statistical mixture of an ensemble of all possible outcomes, which is referred to as non-referring measurement process Ueda . The master equation describing the state evolution in the non-referring measurement is obtained as follows. In the non-referring measurement process, $`\widehat{\rho }(t+dt)`$ is a statistical mixture of the density operator $`\widehat{\rho }_1(t+dt)`$ of the one-count process and the density operator $`\widehat{\rho }_0(t+dt)`$ of the no-count process:

$$\widehat{\rho }(t+dt)=\gamma \widehat{n}_t\widehat{\rho }_1(t+dt)dt+\left(1\gamma \widehat{n}_tdt\right)\widehat{\rho }_0(t+dt).$$
(13)
Substituting $`\widehat{\rho }_1(t+dt)`$ for $`\widehat{\rho }(t_+)`$ in Eq. (5) gives

$$\gamma \widehat{n}_t\widehat{\rho }_1(t+dt)dt=\gamma \widehat{a}\widehat{\rho }(t)\widehat{a}^{}dt.$$
(14)
Setting $`T=dt`$ in Eq. (9) gives

$$\widehat{\rho }_0(t+dt)=\widehat{\rho }(t)+\frac{i}{\mathrm{}}[\widehat{\rho }(t)\widehat{H}_{\mathrm{nc}}^{}\widehat{H}_{\mathrm{nc}}\widehat{\rho }(t)]dt+\gamma \widehat{n}_t\widehat{\rho }(t)dt+O(dt^2).$$
(15)
Substituting Eqs. (14) and (15) into Eq. (13) gives the master equation for the non-referring measurement:

$$\frac{\widehat{\rho }(t)}{t}=\frac{i}{\mathrm{}}[\widehat{\rho }(t),\widehat{H}_0^{\mathrm{rot}}]+\frac{\gamma }{2}\left[2\widehat{a}\widehat{\rho }(t)\widehat{a}^{}\widehat{n}\widehat{\rho }(t)\widehat{\rho }(t)\widehat{n}\right].$$
(16)
When the cavity contains no atoms or when $`g=0`$, we can show that the state of the photon field eventually evolves into a coherent state $`|\alpha `$ with $`\alpha =A/(\delta _{\mathrm{cp}}i\gamma /2)`$ in both referring and non-referring measurement processes, unless the initial state is orthogonal with $`|\alpha ^{}`$ (see Appendix B). This fact is interesting in that the photon number never increases for the steady state $`|\alpha `$ despite the fact that the cavity is being pumped and that no photons are detected. This can be understood from the equation describing the time evolution of the mean photon number in the no-count process

$$\frac{d\widehat{n}_t}{dt}=iA\widehat{a}\widehat{a}^{}_t\gamma (\mathrm{\Delta }\widehat{n})^2_t,$$
(17)
which can be obtained from Eq. (10). For the steady state $`|\alpha `$ with $`\alpha =A/(\delta _{\mathrm{cp}}i\gamma /2)`$, the pumping term \[the first term of the rhs of Eq. (17)\] balances with the last term in Eq. (17) which describes the decrease in the number of photons due to the state reduction by the no-count process.

## III Generation of Squeezed Spin State by Continuous Photodetection

### III.1 Squeezed spin state

Because the spin operators $`\widehat{S}_j`$ ($`j=x,y,z`$) obey the commutation relation $`[\widehat{S}_x,\widehat{S}_y]=i\widehat{S}_x`$ and the corresponding cyclic permutations, the spin vector fluctuates around the mean value $`\widehat{𝐒}=(\widehat{S}_x,\widehat{S}_y,\widehat{S}_z)`$. For example, let us consider the lowest eigenstate of $`\widehat{S}_z`$, $`|S,S`$, where $`S(S+1)`$ is the eigenvalue of $`\widehat{𝐒}^2\widehat{S}_x^2+\widehat{S}_y^2+\widehat{S}_z^2`$. The mean spin vector of this state is given by $`\widehat{𝐒}=(0,0,S)`$. Quantum fluctuations around this are isotropic and given by $`(\mathrm{\Delta }\widehat{S}_x)^2=(\mathrm{\Delta }\widehat{S}_y)^2=S/2`$, satisfying the minimum uncertainty relation $`(\mathrm{\Delta }\widehat{S}_x)^2(\mathrm{\Delta }\widehat{S}_y)^2=\widehat{S}_z^2/4`$. If the fluctuation of one component perpendicular to the mean spin vector, say $`(\mathrm{\Delta }\widehat{S}_x)^2`$, is suppressed below $`|\widehat{S}_z|/2`$, the spin state is defined as squeezed Kitagawa .

A general definition of spin squeezing is obtained by rotating the above situation in spin space, when the mean spin vector $`\widehat{𝐒}`$ points in an arbitrary direction. The definition is thus given by Kitagawa

$$\zeta \frac{2(\mathrm{\Delta }\widehat{S}_{})^2}{|\widehat{𝐒}|}<1,$$
(18)
where $`\widehat{S}_{}`$ is one spin component perpendicular to the mean spin vector, $`\mathrm{\Delta }\widehat{S}_{}\widehat{S}_{}\widehat{S}_{}`$, and $`|\widehat{𝐒}|^2=\widehat{S}_x^2+\widehat{S}_y^2+\widehat{S}_z^2`$. The squeezing parameter $`\zeta `$ is unity in the coherent spin state

$$|\theta ,\varphi =\underset{j}{}\left(\mathrm{cos}\frac{\theta }{2}|\mathrm{e}_j+e^{i\varphi }\mathrm{sin}\frac{\theta }{2}|\mathrm{g}_j\right).$$
(19)
The definition (18) of spin squeezing can also be applied to situations in which the spin state is a mixed state. The physical meaning of the SSS in atomic systems is discussed in Ref. SaitoA .

For an $`N`$-atom system, there are $`N/2+1`$ (when $`N`$ is even) or $`(N+1)/2`$ (when $`N`$ is odd) subspaces that can be indexed by the value of the total spin $`S=N/2,N/21,\mathrm{},0(1/2)`$. In the present paper, we restrict ourselves to the subspace of $`S=N/2`$. This corresponds to the case in which the initial state is the coherent spin state, which is easily prepared in experiments. Since the Hamiltonian (4) commutes with $`\widehat{𝐒}^2`$, the subspaces do not mix by time evolution.

### III.2 Analytic results

For simplicity, we consider the case of $`\delta _{\mathrm{cp}}=\delta _{\mathrm{ap}}=0`$, i.e., the case in which both the coherent pumping field and the atomic transition energy are resonant with the cavity mode. The non-Hermitian operator $`\widehat{H}_{\mathrm{nc}}`$ describing the no-count process (9) is then given by

$`\widehat{H}_{\mathrm{nc}}/\mathrm{}`$ $`=`$ $`A(\widehat{a}+\widehat{a}^{})+g(\widehat{a}\widehat{S}_++\widehat{a}^{}\widehat{S}_{}){\displaystyle \frac{i\gamma }{2}}\widehat{n}`$ (20)
$``$ $`\widehat{a}\widehat{\sigma }^{}+\widehat{a}^{}\widehat{\sigma }{\displaystyle \frac{i\gamma }{2}}\widehat{n},`$
where we define $`\widehat{\sigma }g\widehat{S}_{}+A`$.

As will be shown in Sec. III.3, the generation of SSS is effective for $`\gamma /gS`$, and hence we derive an approximate master equation of the no-count process for large $`\gamma `$. In the following derivation, we assume $`\gamma A`$ and $`AgS`$. Defining $`\widehat{\rho }^{}(t)e^{\gamma t\widehat{n}/2}\widehat{\rho }(t)e^{\gamma t\widehat{n}/2}`$, we rewrite the master equation (10) for the no-count process as

$`{\displaystyle \frac{\widehat{\rho }^{}(t)}{t}}`$ $`=`$ $`i\widehat{\rho }^{}(t)\left(e^{\frac{\gamma t}{2}}\widehat{a}\widehat{\sigma }^{}+e^{\frac{\gamma t}{2}}\widehat{a}^{}\widehat{\sigma }\right)`$
$`i\left(e^{\frac{\gamma t}{2}}\widehat{a}\widehat{\sigma }^{}+e^{\frac{\gamma t}{2}}\widehat{a}^{}\widehat{\sigma }\right)\widehat{\rho }^{}(t)+\gamma \widehat{n}_t\widehat{\rho }^{}(t).`$
Since $`\gamma A`$, the photon field is almost always in the vacuum state, which is justified below. We may therefore approximate $`\widehat{\rho }^{}(t)`$ as $`\widehat{\rho }^{}(t)|00|\widehat{\rho }_s(t)`$, where $`\widehat{\rho }_s(t)`$ denotes the reduced density operator for the spin. Integrating Eq. (III.2) iteratively gives

$`\widehat{\rho }^{}(t+\mathrm{\Delta }t)`$ $`=`$ $`|00|\{\widehat{\rho }_s(t)+\gamma {\displaystyle _t^{t+\mathrm{\Delta }t}}dt_1\widehat{n}_{t_1}\widehat{\rho }_s(t)`$ (22)
$`{\displaystyle _t^{t+\mathrm{\Delta }t}}dt_2{\displaystyle _t^{t_2}}dt_1[e^{\frac{\gamma (t_1t_2)}{2}}(\widehat{\sigma }^{}\widehat{\sigma }\widehat{\rho }_s(t)+\widehat{\rho }_s(t)\widehat{\sigma }^{}\widehat{\sigma })+\gamma ^2\widehat{n}_{t_1}\widehat{n}_{t_2}]\}`$
$`+|11|{\displaystyle _t^{t+\mathrm{\Delta }t}}𝑑t_2{\displaystyle _t^{t_2}}𝑑t_12e^{\frac{\gamma (t_1+t_2)}{2}}\widehat{\sigma }\widehat{\rho }_s(t)\widehat{\sigma }^{}+\mathrm{},`$
where we have kept only relevant terms. We take $`\mathrm{\Delta }t`$ such that $`\gamma \mathrm{\Delta }t1`$ and $`\gamma \mathrm{\Delta }tA^2/\gamma ^21`$ are satisfied. From the last term of Eq. (22), we find that the mean photon number is given by $`\widehat{n}_{t+\mathrm{\Delta }t}4\gamma ^2\widehat{\sigma }^{}\widehat{\sigma }_t`$, which is $`O(A^2/\gamma ^2)`$. The last term and the term including $`\gamma ^2\widehat{n}_{t_1}\widehat{n}_{t_2}`$ in Eq. (22) are therefore $`O(A^2/\gamma ^2)`$ and can be ignored. Thus, Eq. (22) can be approximated as

$$\widehat{\rho }(t+\mathrm{\Delta }t)|00|\left\{\widehat{\rho }_s(t)\frac{2\mathrm{\Delta }t}{\gamma }[\widehat{\sigma }^{}\widehat{\sigma }\widehat{\rho }_s(t)+\widehat{\rho }_s(t)\widehat{\sigma }^{}\widehat{\sigma }]+\frac{4\mathrm{\Delta }t}{\gamma }\widehat{\sigma }^{}\widehat{\sigma }_t\widehat{\rho }_s(t)\right\}.$$
(23)
Since, the last two terms in the curly brackets in Eq. (23) are $`O(\gamma \mathrm{\Delta }tA^2/\gamma ^2)`$, which is assumed to be small, we obtain an effective master equation of the spin state for the no-count process:

$$\frac{\widehat{\rho }_s(t)}{t}=\frac{2}{\gamma }[\widehat{\sigma }^{}\widehat{\sigma }\widehat{\rho }_s(t)+\widehat{\rho }_s(t)\widehat{\sigma }^{}\widehat{\sigma }]+\frac{4}{\gamma }\widehat{\sigma }^{}\widehat{\sigma }_t\widehat{\rho }_s(t).$$
(24)
It follows from Eq. (24) that the system evolves into the lowest eigenstate of the operator

$$\widehat{\sigma }^{}\widehat{\sigma }=g^2\widehat{S}_+\widehat{S}_{}+gA(\widehat{S}_++\widehat{S}_{})+A^2,$$
(25)
if the initial spin state is not orthogonal with the lowest eigenstate. To see this, let us expand the density operator as

$$\widehat{\rho }(t)=\underset{n,m}{}c_{nm}|\varphi _n\varphi _m|,$$
(26)
where $`|\varphi _n`$ is the eigenstate of the operator (25) with eigenvalue $`\epsilon _n`$, i.e., $`\widehat{\sigma }^{}\widehat{\sigma }|\varphi _n=\epsilon _n|\varphi _n`$. The master equation (24) is rewritten as

$$\dot{c}_{nm}(t)=\frac{2}{\gamma }(2\widehat{\sigma }^{}\widehat{\sigma }_t\epsilon _n\epsilon _m)c_{nm}(t),$$
(27)
which can be solved to give

$$c_{nm}(t)=\mathrm{exp}\left[\frac{4}{\gamma }^t𝑑t^{}\widehat{\sigma }^{}\widehat{\sigma }_t^{}\frac{2}{\gamma }(\epsilon _n+\epsilon _m)t\right]c_{nm}(0).$$
(28)
Thus, the system develops into the state that minimizes $`\epsilon _n+\epsilon _m`$, that is, the lowest eigenstate of the operator (25), if the corresponding initial coefficient $`c_{nm}(0)`$ is not zero. We note that Eq. (25) is independent of $`\gamma `$, and therefore small experimental fluctuations in $`\gamma `$ do not affect the steady state. Since the relaxation coefficient in Eq. (24) is proportional to $`\gamma ^1`$, the relaxation becomes slower for larger $`\gamma `$. This is because the number of photons in the cavity is smaller for larger $`\gamma `$, and therefore the state reduction caused by the no-count process becomes less effective.

The steady state of the master equation (24) is the SSS, as will be numerically shown in Sec. III.3. We evaluate perturbatively the degree of spin squeezing (18) of the steady state. When $`Ag`$, the first-order perturbation with respect to $`g/A`$ yields (Appendix C)

$$\zeta =1\frac{g}{4A}(N1)<1,$$
(29)
showing that the state is an SSS. When $`Ag`$, the second-order perturbation with respect to $`A/g`$ yields (Appendix C)

$$\zeta =1\frac{2A^2}{N^2g^2}<1,$$
(30)
showing that the state is again an SSS.

### III.3 Exact numerical results

We numerically simulate the stochastic process of the continuous measurement described by Eqs. (5) and (10) with $`\widehat{H}_{\mathrm{nc}}`$ given by Eq. (20). We divide the time into a small interval $`\mathrm{\Delta }t`$ such that $`\gamma \widehat{n}_t\mathrm{\Delta }t`$ is always much smaller than unity. We then generate a random number $`r`$ distributed uniformly between 0 and 1 for each time interval. If $`0r\gamma \widehat{n}_t\mathrm{\Delta }t`$, we assume that the one-count process occurs between $`t`$ and $`t+\mathrm{\Delta }t`$ and calculate the density operator at $`t+\mathrm{\Delta }t`$ according to Eq. (5). If $`\gamma \widehat{n}_t\mathrm{\Delta }t<r1`$, we assume that the no-count process occurs during $`t`$ and $`t+\mathrm{\Delta }t`$ and calculate the density operator at $`t+\mathrm{\Delta }t`$ according to Eq. (10).

Figure 2 shows an example of the stochastic time evolutions of various quantities in the referring measurement process with $`A/g=5`$, $`\gamma /g=10`$, and $`N=10`$, where the photon field is initially in the vacuum state and the atoms are in $`|S,S`$.

The impulses that can be seen at the bottom of Fig. 2 (c) indicate the times at which photons are detected. They are clearly bunched, reflecting the fact that $`g^{(2)}\widehat{a}^2\widehat{a}^2/\widehat{n}^2`$ in the cavity is larger than unity most of the time as shown in Fig. 2 (a). The bunching nature of photons makes the duration of each no-count process longer, helping the system reach the no-count steady state (shown as plateaus in Fig. 2). We note that in the no-count steady state, the SSS is generated with $`\zeta 0.56`$ \[solid curve in Fig. 2 (c)\], where $`\zeta `$ is defined in Eq. (18). In the non-referring measurement process, in contrast, no SSS is obtained \[dashed line in Fig. 2 (c), obtained by numerically integrating Eq. (16)\], indicating that the information concerning the detection times of photons is crucial for obtaining the SSS. It should be emphasized that in such a referring measurement process $`\zeta `$ is smaller than unity most of the time, and therefore the SSS is obtained as a quasi-steady state.

Figure 3 (a) shows the quasi-probability distribution (QPD) corresponding to the spin state at $`t=15`$, when the no-count steady state is reached.

We can see that the spin fluctuation is squeezed in the azimuthal direction. However, at $`t=16`$ the spin state is disturbed by the photon counts and the corresponding QPD shown in Fig. 3 (b) exhibits a pattern similar to a rotated Dicke state.

The degree of spin squeezing, $`\zeta `$, depends on the parameters $`A`$, $`\gamma `$, and $`N`$. Figure 4 shows the dependence of $`\zeta _{\mathrm{st}}`$ on $`A`$ and $`\gamma `$ for $`N=2`$ and $`10`$, where $`\zeta _{\mathrm{st}}`$ is the steady-state value of $`\zeta `$ in the no-count process.

We see that the SSS is obtained for a wide range of parameters. We find that in both cases of $`N=2`$ and $`N=10`$ the maximum squeezing is obtained for $`A/gN/2=S`$, and we have numerically confirmed that this finding holds true for other values of $`N`$. For $`\gamma A`$, $`\zeta _{\mathrm{st}}`$ depends little on $`\gamma `$, in agreement with the analytic results of Eqs. (25), (29), and (30).

Figure 5 shows the dependence of $`\zeta _{\mathrm{st}}`$ and the angle of the spin vector $`\theta _{\mathrm{st}}`$ on $`N`$, where the squares are obtained by solving the master equation (10) and the circles by diagonalizing Eq. (25).

Both plots are in excellent agreement, confirming the validity of the analytical treatment in Sec. III.2. It is interesting to note that both $`\zeta _{\mathrm{st}}`$ and $`\theta _{\mathrm{st}}`$ are proportional to $`N^{1/3}`$. The other models of spin squeezing also exhibit such power laws: $`N^{2/3}`$ for one-axis twisting and $`N^1`$ for two-axis counter twisting Kitagawa . In a high-Q cavity, the Jaynes-Cummings interaction between photons and atoms generates spin squeezing proportional to $`N^{1/4}`$ SaitoA , demonstrating that the present scheme is more effective for larger $`N`$. Table I summarizes the dependence of $`\zeta `$ on $`N`$ in these schemes.

## IV Discussions and Conclusions

When the number of excited atoms is much smaller than the total number of atoms $`N`$, the Holstein-Primakoff transformation Holstein  can be approximated by $`\widehat{S}_{}(2S)^{1/2}\widehat{b}`$ and $`\widehat{S}_+(2S)^{1/2}\widehat{b}^{}`$, where $`\widehat{b}`$ and $`\widehat{b}^{}`$ are the bose operators satisfying $`[\widehat{b},\widehat{b}^{}]=1`$. The non-Hermitian operator describing the no-count process (20) then becomes

$$\widehat{H}_{\mathrm{nc}}/\mathrm{}=A(\widehat{a}+\widehat{a}^{})+G(\widehat{a}\widehat{b}^{}+\widehat{a}^{}\widehat{b})\frac{i\gamma }{2}\widehat{n},$$
(31)
where $`G(2S)^{1/2}g`$. In a similar manner to the derivation given in Appendix B, one can show that during the no-count process the system develops into $`|0_a|\beta _b`$, where $`|0_a`$ is the vacuum state of the photon field and $`|\beta _b=e^{\beta (\widehat{b}^{}\widehat{b})}|0_b`$ is the coherent state of the boson with $`\beta =A/G`$. Therefore, the bosonization of spin eliminates the squeezing properties of atoms, and thus the spin algebra, which plays a role in nonlinear saturation, is crucial for obtaining squeezing in the present scheme.

Experimentally, the model presented may be implemented in various systems, such as a microwave cavity with Rydberg atoms Haroche , an optical cavity with cold atoms Kimble , and trapped ions coupled with their center-of-mass motion Wineland . For example, let us consider a system in which Rb atoms in the circular Rydberg state pass through a superconducting microwave cavity Haroche . The coupling constant between the atoms and the cavity photons can be varied from 0 to 24 kHz Haroche . If we take $`\kappa /(2\pi )=2`$ kHz, $`g/(2\pi )=100`$ Hz, $`\tau _\mathrm{p}=10`$ $`\mu \mathrm{s}`$, and $`N=10`$, the conditions (34) are satisfied as $`\kappa \tau _\mathrm{p}0.1`$ and $`gS\tau _\mathrm{p}A\tau _\mathrm{p}0.03`$, where $`\kappa `$ is the coupling constant between probe atoms and photons, and $`\tau _\mathrm{p}`$ is the passage time of a probe atom through the cavity (see Appendix A). These parameters give $`\gamma /g=\kappa ^2\tau _\mathrm{p}/g2.5`$, and the lifetime of the Rydberg state $`30`$ ms Haroche  imposes the constraint $`gt<19`$. The cavity relaxation time now reaches 0.3 s Walther . Thus the scheme is experimentally feasible.

In conclusion, we have studied the system in which $`N`$ two-level atoms interact with photons in a pumped cavity, and have shown that the atomic state develops into an SSS by continuous measurement of photon number in the cavity. In particular, the SSS is obtained as a quasi-steady state during the no-count process. The present method does not require fine tuning of the initial state and the time and strength of interaction, and therefore significantly facilitates experimental realization of the SSS.

###### Acknowledgements.

This work was supported by the Special Coordination Funds for Promoting Science and Technology and a Grant-in-Aid for Scientific Research (Grant No. 11216204) by the Ministry of Education, Science, Sports, and Culture of Japan, by the Yamada Science Foundation, and by the Toray Science Foundation.

## Appendix A A model of the continuous measurement of photon number

For continuous measurement of photon number, the Srinivas-Davies “ansatz” Srinivas  is usually employed. However, it is by no means obvious whether it is valid in our system, and therefore, we explicitly derive nonunitary time evolutions (5) and (9) using the continuous photodetection model of Ref. Imoto .

We consider a situation in which probe atoms prepared in the ground state $`|\mathrm{g}`$ pass through the cavity successively, and the internal states of the atoms coming out of the cavity are detected by, e.g., an ionization detector (Fig. 1). When a probe atom is detected to be in the excited state $`|\mathrm{e}`$, we regard this event as a one-count process. When the atom coming out of the cavity remains in the ground state, we regard this event as a no-count process.

The time evolution as the probe atom passes through the cavity is described by the Hamiltonian

$$\widehat{H}_\mathrm{p}=\widehat{H}_0^{\mathrm{rot}}+\mathrm{}\kappa (\widehat{a}\widehat{\sigma }_++\widehat{a}^{}\widehat{\sigma }_{}),$$
(32)
where $`\widehat{H}_0^{\mathrm{rot}}`$ is given in Eq. (4), $`\widehat{\sigma }_+`$ and $`\widehat{\sigma }_{}`$ are the raising and lowering operators of the probe atom, and $`\kappa `$ is a coupling constant which we assume to be real. We assume that the transition frequency of the probe atom is resonant with the cavity mode. When the probe atom enters the cavity at time $`t`$ and comes out of it at $`t+\tau _\mathrm{p}`$, the density operator $`\widehat{\rho }_{\mathrm{sp}}`$ for the system and the probe atom evolves as

$`\widehat{\rho }_{\mathrm{sp}}(t+\tau _\mathrm{p})`$ $`=`$ $`\widehat{\rho }(t)|\mathrm{g}\mathrm{g}|{\displaystyle \frac{i\tau _\mathrm{p}}{\mathrm{}}}[\widehat{H}_\mathrm{p},\widehat{\rho }(t)|\mathrm{g}\mathrm{g}|]`$ (33)
$`{\displaystyle \frac{\tau _\mathrm{p}^2}{2\mathrm{}^2}}[\widehat{H}_\mathrm{p}[\widehat{H}_\mathrm{p},\widehat{\rho }(t)|\mathrm{g}\mathrm{g}|]]+\mathrm{},`$
where $`\widehat{\rho }(t)`$ is the density operator for the system (excluding the probe atom).

Here we assume that $`\tau _\mathrm{p}`$ is sufficiently small so that third and higher-order terms in Eq. (33) can be ignored. This approximation is valid if the change in the state of the system is small during $`\tau _\mathrm{p}`$. This condition is satisfied if the following conditions are met:

$$A\tau _\mathrm{p}^11,gS\tau _\mathrm{p}^11,\delta _{\mathrm{cp}}\tau _\mathrm{p}^11,\delta _{\mathrm{ap}}\tau _\mathrm{p}^11.$$
(34)
When the probe atom is detected to be in the excited state $`|\mathrm{e}`$, the density operator of the system is projected into

$$\mathrm{Tr}_{\mathrm{probe}}[|\mathrm{e}\mathrm{e}|\widehat{\rho }_{\mathrm{sp}}(t+\tau _\mathrm{p})]=\kappa ^2\tau _\mathrm{p}^2\widehat{a}\widehat{\rho }(t)\widehat{a}^{}.$$
(35)
The normalized form of Eq. (35) gives the density operator immediately after the no-count process:

$$\widehat{\rho }(t_+)=\frac{\widehat{a}\widehat{\rho }(t)\widehat{a}^{}}{\mathrm{Tr}[\widehat{a}\widehat{\rho }(t)\widehat{a}^{}]},$$
(36)
where $`t_+`$ denotes a time infinitesimally later than $`t`$. The probability that the one-count process (35) occurs is given by the trace of Eq. (35), i.e., $`\gamma \tau _\mathrm{p}\widehat{n}_t`$, where

$$\gamma \kappa ^2\tau _\mathrm{p}$$
(37)
can be interpreted as the rate of photodetection per one photon.

When the probe atom is detected in its ground state $`|g`$, the density operator of the system is projected into

$`\mathrm{Tr}_{\mathrm{probe}}[|\mathrm{g}\mathrm{g}|\widehat{\rho }_{\mathrm{sp}}(t+\tau _\mathrm{p})]`$ (38)
$``$ $`\widehat{\rho }(t){\displaystyle \frac{i\tau _\mathrm{p}}{\mathrm{}}}[\widehat{H}_0^{\mathrm{rot}},\widehat{\rho }(t)]{\displaystyle \frac{\tau _\mathrm{p}^2}{2\mathrm{}^2}}[\widehat{H}_0^{\mathrm{rot}},[\widehat{H}_0^{\mathrm{rot}},\widehat{\rho }(t)]`$
$`{\displaystyle \frac{\tau _\mathrm{p}^2\kappa ^2}{2}}\left[\widehat{n}\widehat{\rho }(t)+\widehat{\rho }(t)\widehat{n}\right]`$
$``$ $`e^{\tau _\mathrm{p}\left(\frac{i}{\mathrm{}}\widehat{H}_0^{\mathrm{rot}}\frac{\gamma }{2}\right)}\widehat{\rho }(t)e^{\tau _\mathrm{p}\left(\frac{i}{\mathrm{}}\widehat{H}_0^{\mathrm{rot}}\frac{\gamma }{2}\right)}.`$
When successive probe atoms are found to be in the ground state over a time interval $`T`$, the density operator of the system after such a no-count process is given by Eq. (9).

Probe atoms were introduced to explicitly derive Eqs. (5) and (9). However, these expressions may also describe a more general system under continuous photodetection, such as the situation in which photons leaking out of the cavity are detected by a photodetector.

## Appendix B Steady state in a pumped cavity under continuous photodetection

We show that the state of the system described by the Hamiltonian (4) eventually evolves into a coherent state in both referring and non-referring measurements.

The non-Hermitian operator that describes the no-count process (9) is rewritten in the rotating frame by

$`\widehat{H}_{\mathrm{nc}}/\mathrm{}`$ $`=`$ $`A(\widehat{a}+\widehat{a}^{})+(\delta _{\mathrm{cp}}i\gamma /2)\widehat{n}`$ (39)
$`=`$ $`\mathrm{\Omega }\left(\widehat{a}^{}+{\displaystyle \frac{A}{\mathrm{\Omega }}}\right)\left(\widehat{a}+{\displaystyle \frac{A}{\mathrm{\Omega }}}\right){\displaystyle \frac{A^2}{\mathrm{\Omega }}}`$
$`=`$ $`\mathrm{\Omega }e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}\widehat{n}e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}{\displaystyle \frac{A^2}{\mathrm{\Omega }}},`$
where we define $`\mathrm{\Omega }\delta _{\mathrm{cp}}i\gamma /2`$. We then obtain

$`e^{i\widehat{H}_{\mathrm{nc}}t/\mathrm{}}`$ $``$ $`e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}e^{i\mathrm{\Omega }t\widehat{n}}e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}`$
$`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}e^{in\delta _{\mathrm{cp}}tn\gamma t/2}e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}|nn|e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}.`$
For $`\gamma t1`$, the term corresponding to $`n=0`$, i.e., $`e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}|00|e^{\frac{A}{\mathrm{\Omega }}(\widehat{a}+\widehat{a}^{})}|\alpha \alpha ^{}|`$ becomes dominant in the summation, where $`|\alpha `$ is the coherent state with $`\alpha =A/\mathrm{\Omega }`$. We thus obtain

$$\underset{t\mathrm{}}{lim}\frac{e^{i\widehat{H}_{\mathrm{nc}}t/\mathrm{}}|\psi }{e^{i\widehat{H}_{\mathrm{nc}}t/\mathrm{}}|\psi }=|\alpha $$
(41)
for $`\alpha ^{}|\psi 0`$. Therefore, during the no-count process (9) the state develops into $`|\alpha `$ unless the initial state is orthogonal to $`|\alpha ^{}`$. If $`\alpha ^{}|\psi =0`$, the term of $`n=1`$ in the summation of Eq. (B) becomes dominant, and the state develops into a displaced number state.

If the photons are detected at times $`t_1,t_2,\mathrm{},t_n`$, the time development operator (12) for the referring measurement takes the form

$`\widehat{P}(T;t_1,\mathrm{},t_n)`$ (42)
$`=`$ $`e^{i\widehat{H}_{\mathrm{nc}}T/\mathrm{}}\mathrm{\Pi }_{k=1}^n\left[e^{i\mathrm{\Omega }t_k}\widehat{a}+{\displaystyle \frac{A}{\mathrm{\Omega }}}(e^{i\mathrm{\Omega }t_k}1)\right].`$
From Eqs. (41) and (42) we obtain

$$\underset{T\mathrm{}}{lim}\widehat{P}(T;t_1,\mathrm{},t_n)|\psi |\alpha .$$
(43)
Thus, the state develops to $`|\alpha `$ in the referring measurement. This result is also true for non-referring measurement, since Eq. (43) depends neither on the number of detected photons nor on the times at which photons are detected.

## Appendix C Proof of Eqs. (29) and (30)

We calculate perturbatively Sakurai  the squeezing parameter for the lowest eigenstate of the operator (25).

When $`g/A1`$, the first term on the rhs of Eq. (25) can be treated as a perturbation. For convenience, we rotate in the spin space by $`\pi /2`$ with respect to the $`S_y`$ axis. The operator (25) then becomes

$$\widehat{\sigma }^{}\widehat{\sigma }=A^2\left[1+\frac{2g}{A}\widehat{S}_z+\frac{g^2}{A^2}(\widehat{𝐒}^2\widehat{S}_x^2\widehat{S}_x)\right].$$
(44)
The unperturbed ground state is then taken to be the eigenstate of $`\widehat{S}_z`$ with eigenvalue $`S`$, which we write as $`|S,S`$. Performing first-order perturbation with respect to the last term of Eq. (44), we obtain the ground state as

$$|S,S+\frac{\sqrt{2S}g}{4A}\left(|S,1S+\frac{1}{4}\sqrt{2(2S1)}|S,2S\right),$$
(45)
with the squeezing parameter calculated to be

$$\zeta =1\frac{g}{4A}(2S1).$$
(46)
When $`Ag`$, the second term on the rhs of Eq. (25) may be treated as a perturbation. In this case, we must perform second-order perturbation to obtain spin squeezing. The ground state is then given by

$`\left(1{\displaystyle \frac{A^2}{4Sg^2}}\right)|S,S{\displaystyle \frac{A}{\sqrt{2S}g}}|S,1S`$
$`+{\displaystyle \frac{A^2}{g^2}}{\displaystyle \frac{1}{\sqrt{4S(2S1)}}}|S,2S,`$ (47)
with the squeezing parameter calculated to be

$$\zeta =1\frac{A^2}{2S^2g^2}.$$
(48)


format: text

text:
# P- and T-violating Schiff moment of the Mercury nucleus

## 1 Introduction

The most precise limit on parity- and time-invariance violating nucleon-nucleon interaction has been obtained from the measurement of the atomic electric dipole moment of 199Hg . The hadronic part of the atomic dipole moment associated with the electric dipole moment of the 199Hg nucleus manifests itself through the Schiff moment which is the first nonzero term in the expansion of the nuclear electromagnetic potential after including the screening of the atomic electrons .

The operator for the Schiff moment is

$$S_\mu =\frac{1}{10}\sqrt{\frac{4\pi }{3}}\underset{i}{\overset{A}{}}e_i\left(r_i^3\frac{5}{3}r^2_{ch}r_i\right)Y_{1\mu }(\widehat{𝐫}_i),$$
(1)
where $`e_i`$ is $`e`$ for a proton and zero for a neutron. The Schiff moment generates T- and P-odd electrostatic potential in the form

$$\varphi (𝐫)=4\pi 𝐒\mathbf{}\delta (𝐫).$$
(2)
Interaction of atomic electrons with the potential given by Eq. (2) produces an atomic dipole moment

$$d_{atom}=\underset{n}{}\frac{0|e_i^Z\varphi (𝐫_i)|nn|e_i^Zz_i|0}{E_nE_0}+\mathrm{h}.\mathrm{c}.$$
(3)
Due to contact origin of the potential, the electrons in $`s`$\- and $`p`$\- atomic orbitals only contribute to the dipole moment given by Eq. (3).

The nuclear Schiff moment was calculated so far in a simplified model  without considering many-body nuclear structure effects. These effects have to be understood properly if we intend to extract the parameters of P- and T-violating nuclear interaction from the value of the Schiff moment. For light nuclei the properties of the Schiff moment strength obtained in modern shell model calculations were discussed in Ref. . The study of the polarization effects associated with the coupling to isoscalar dipole compression mode has been performed in Ref. . In our paper we calculate the Schiff moment of 199Hg nucleus within RPA framework with effective residual strong forces using finite range P- and T-odd weak nuclear interaction.

## 2 Basic ingredients of the theory

### 2.1 Nucleon-nucleon P- and T-odd interaction

We use the interaction generated by P- and T-violating pion exchange .

$$W(𝐫_1𝐫_2)=\frac{g}{8\pi m_p}[(g_0𝝉_1𝝉_2+g_2(𝝉_1𝝉_23\tau _1^3\tau _2^3))(𝝈_1𝝈_2)$$

$$+g_1(\tau _1^3𝝈_1\tau _2^3𝝈_2)]\mathbf{}_1\frac{e^{m_\pi r_{12}}}{r_{12}},$$
(4)
where $`g`$ is the usual strong pion-nucleon pseudoscalar coupling constant, $`g_0`$, $`g_1`$, and $`g_2`$ correspond to isoscalar, isovector , and isotensor P- and T-odd couplings, $`m_p`$ is the proton mass. In contrast to P-odd and T-even interaction, in Eq. (4) the exchange of $`\pi ^0`$ is allowed. This term produces the direct contribution to P- and T- odd part of nuclear mean field while the other terms produce the exchange contribution only. Since the direct contribution dominates for finite range potentials we can expect that the interaction (4) is the leading one and the exchange of heavier mesons can be omitted.

In previous calculations the phenomenological contact interaction has often been used instead of finite range interaction given by Eq. (4). It has the form

$$W_c(𝐫_a𝐫_b)=\frac{G}{\sqrt{2}}\frac{1}{2m_p}((\eta _{ab}𝝈_a\eta _{ba}𝝈_b)\text{ }\mathbf{}_a\delta (𝐫_a𝐫_b)+$$

$$\eta _{ab}^{}[𝝈_a\times 𝝈_b]\{(𝐩_a𝐩_b),\delta (𝐫_a𝐫_b)\}),$$
(5)
where $`G`$ is the Fermi constant. In the limit $`m_\pi \mathrm{}`$ the interaction (4) transforms into Eq. (5) after the substitution $`gg_i\frac{Gm_\pi ^2}{\sqrt{2}}\eta `$. We shall use this factor when comparing our results with those obtained using the contact interaction given by Eq. (5).

### 2.2 Nuclear mean field and correction from the weak forces

In our calculations we used full single-particle spectrum including continuum. The single-particle basis was obtained using partially self-consistent mean-field potential of Ref. . The potential includes four terms. The isoscalar term is the standard Woods-Saxon potential

$$U_0(r)=\frac{V}{1+\mathrm{exp}\frac{rR}{a}},$$
(6)
with the parameters $`V=52.03`$ MeV, $`R=1.2709A^{1/3}`$ fm, and $`a=0.742`$ fm. Two other terms $`U_{ls}(r)`$, and $`U_\tau (r)`$ were calculated self-consistently using two-body Migdal-type interaction of Ref.  for the spin-orbit and isovector parts of the potential. The last term is the Coulomb potential calculated for uniformly charged sphere with $`R_C=1.18A^{1/3}`$ fm. The mean field potential obtained in this way produces good fit for single particle energies and r.m.s. radii for nuclei in the lead region.

The correction to the mean field (6) from the weak interaction (4) consist of direct and exchange terms. The direct term has the form

$$\delta U_{dir}(𝐫)=\frac{gm_\pi ^2}{\pi m_p}(𝝈𝐧)\tau ^3_0^{\mathrm{}}r^2dr^{}b_{\mathrm{1\hspace{0.17em}0}}(r,r^{})[(g_02g_2)(\rho _p(r)\rho _n(r))+$$

$$g_1(\rho _p(r)+\rho _n(r))],$$
(7)
where the function $`b_{\mathrm{1\hspace{0.17em}0}}(r_1,r_2)`$ is a combination of spherical Bessel functions of imaginary argument

$$b_{l_1l_2}(r_1,r_2)=i_{l_1}(m_\pi r_1)k_{l_2}(m_\pi r_2)\theta (r_2r_1)i_{l_2}(m_\pi r_2)k_{l_1}(m_\pi r_1)\theta (r_1r_2).$$
(8)
Note, that the potential given by Eq. (7) is pure isovector. The contribution of the isovector interaction component dominates in Eq. (7), the isoscalar and isotensor components of the interaction are suppressed by the factor $`\frac{NZ}{A}`$. For zero range interaction the potential would be proportional to $`\rho (r)`$. The gradient makes this potential very sensitive to the details of nuclear surface. Our potential given by Eq. (7) is less sensitive to the surface due to additional integration over region of the order of pion Compton wavelength.

The exchange term is more complicated. The matrix element of it taken over angular variable is a non-local operator in radial coordinates.

$$\stackrel{~}{\nu }|\delta U_{exch}(𝐫,𝐫^{})|\nu =W_\nu (r,r^{}),$$

where $`|\stackrel{~}{\nu }=(𝝈𝐧)|\nu `$, and

$$W_\nu (r,r^{})=\frac{1}{2j_\nu +1}\frac{gm_\pi ^2}{\pi m_p}Tr_2\{\underset{\kappa l_1l_2}{}(\frac{g_0}{2}(3\tau _1^3\tau _2^3)2g_2\tau _1^3\tau _2^3+\frac{g_1}{2}(\tau _1^3+\tau _2^3))\times $$

$$n_\kappa \left(\begin{array}{ccc}l_1& l_2& 1\\ 0& 0& 0\end{array}\right)b_{l_1l_2}(r,r^{})R_\kappa (r)R_\kappa (r^{})[()^{l_1}[l_1](\kappa ||T_{l_2}^{l_1}||\stackrel{~}{\nu })^{}(\kappa ||Y_{l_2}||\nu )$$

$$()^{l_2}[l_2](\stackrel{~}{\nu }||Y_{l_1}||\kappa )(\nu ||T_{l_1}^{l_2}||\kappa )^{}]\}.$$
(9)
The Trace is assumed over isospin variable of the second particle , $`n_\kappa `$ is the occupation number of the single particle state $`|\kappa `$, $`[l]=\sqrt{2l+1}`$ , and the tensor operator $`T_{JM}^L(𝐧)=\{𝝈Y_L(𝐧)\}_{JM}`$.

The correction to the single particle wave function $`\psi _\nu (𝐫)`$ can be presented as

$$\delta \psi _\nu (𝐫)=(𝝈𝐧)\mathrm{\Omega }_\nu (𝐧)\delta R_\nu (r),$$

where $`\mathrm{\Omega }_\nu (𝐧)`$ is the angular part of the wave function. The radial correction $`\delta R_\nu (r)`$ is the sum of the direct and the exchange terms $`\delta R_\nu (r)=\delta R_{\nu dir}(r)+\delta R_{\nu exch}(r)`$, where

$$\delta R_{\nu dir}(r)=_0^{\mathrm{}}\text{G}_{j_\nu \stackrel{~}{l}_\nu }(r,r^{}|ϵ_\nu )\delta U_{dir}(r^{})R_\nu (r^{})r^2𝑑r^{},$$

$$\delta R_{\nu exch}(r)=_0^{\mathrm{}}\text{G}_{j_\nu \stackrel{~}{l}_\nu }(r,r^{}|ϵ_\nu )W_\nu (r^{},r^{\prime \prime })R_\nu (r^{\prime \prime })r^2r^{\prime \prime 2}𝑑r^{}𝑑r^{\prime \prime }.$$
(10)
Here $`\text{G}_{j_\nu \stackrel{~}{l}_\nu }(r,r^{}|ϵ_\nu )`$ is the Green function of radial Schrödinger equation for the total angular momentum $`j_\nu `$, and the orbital angular momentum $`\stackrel{~}{l}_\nu =2j_\nu l_\nu `$ . $`ϵ_\nu `$ is the single particle energy.

## 3 Core polarization

The effects of the core polarization for a one particle operator can be treated introducing a renormalized operator $`\stackrel{~}{𝐒}`$ satisfying the equation

$$\stackrel{~}{𝐒}_{\nu ^{}\nu }=𝐒_{\nu ^{}\nu }^0+\underset{\mu ^{}\mu }{}\stackrel{~}{𝐒}_{\mu \mu ^{}}\frac{n_\mu n_\mu ^{}}{ϵ_\mu ϵ_\mu ^{}+\omega }\nu ^{}\mu ^{}|F+W|\mu \nu ,$$
(11)
where $`𝐒^0`$ is the bare Schiff moment operator given by Eq. (1), $`n_\mu `$ and $`ϵ_\mu `$ are the single particle occupation numbers and energies. For static moments the external frequency $`\omega 0`$. The interaction in Eq. (11) includes both the strong residual interaction $`F`$ and the weak one. The latter we take in the form given by Eq. (4). Strictly speaking, the interaction of two nucleons in nuclear matter differs from the interaction in the vacuum. We do not discuss this effect here and keep the weak interaction in the form (4). The single particle wave functions in Eq. (11) are the eigenstates of the mean field which is also the sum of the strong and the weak fields given by Eqs. (7,9). Since the weak forces are really small compared to the strong interaction, it is natural to treat them perturbatively. The simplest way to do it is to present $`<𝐫|\nu >=\psi _\nu (𝐫)+\delta \psi _\nu (𝐫)`$ and to gather the terms linear in $`\delta \psi `$ and $`W(𝐫_1𝐫_2)`$. It is convenient to present the matrix element of the Schiff moment as a sum of three terms

$$\stackrel{~}{𝐒}_{\nu \nu }=\delta \psi _\nu |𝐒|\psi _\nu +\psi _\nu |𝐒|\delta \psi _\nu +\psi _\nu |\delta 𝐒|\psi _\nu .$$
(12)
The operator $`𝐒`$ in the first two terms satisfies the same Eq. (11) where the strong interaction only is included. The corrections from the weak forces for these terms are in the wave functions of an odd nucleon only. The third term presents an induced contribution  arising from the P- and T-violating corrections to the intermediate states $`|\mu `$ and $`|\mu ^{}`$. The equation for $`\delta 𝐒`$ is

$$(\delta 𝐒\delta 𝐒_{NL})_{\nu ^{}\nu }=(\delta 𝐒_0)_{\nu ^{}\nu }+\underset{\mu ^{}\mu }{}(\delta 𝐒\delta 𝐒_{NL})_{\mu \mu ^{}}\frac{n_\mu n_\mu ^{}}{ϵ_\mu ϵ_\mu ^{}}\nu ^{}\mu ^{}|F|\mu \nu ,$$
(13)
where $`\delta 𝐒_{NL}`$ is the non-local part of $`\delta 𝐒`$ produced by the exchange matrix elements of weak interaction $`W`$.

$$(\delta 𝐒_{\mathrm{𝐍𝐋}})_{\nu ^{}\nu }=\underset{\mu \mu ^{}}{}(𝐒)_{\mu \mu ^{}}\frac{n_\mu n_\mu ^{}}{ϵ_\mu ϵ_\mu ^{}}\nu ^{}\mu ^{}|W_{exch}|\mu \nu .$$
(14)
The equation for $`\delta 𝐒_0`$ is

$$(\delta 𝐒_0)_{\nu ^{}\nu }=\underset{\mu \mu ^{}}{}\frac{n_\mu n_\mu ^{}}{ϵ_\mu ϵ_\mu ^{}}(\mu |𝐒|\mu ^{}\nu ^{}\mu ^{}|W_{dir}|\mu \nu +\mu |\delta 𝐒_{NL}|\mu ^{}\nu ^{}\mu ^{}|F|\mu \nu +$$

$$\delta \psi _\mu |𝐒|\mu ^{}\nu ^{}\mu ^{}|F|\mu \nu +\mu |𝐒|\delta \psi _\mu ^{}\nu ^{}\mu ^{}|F|\mu \nu +$$

$$\mu |𝐒|\mu ^{}\nu ^{}\delta \psi _\mu ^{}|F|\mu \nu +\mu |𝐒|\mu ^{}\nu ^{}\mu ^{}|F|\delta \psi _\mu \nu ).$$
(15)
Note, that in the absence of the core polarization the first term only in Eq. (15) contributes into $`\delta 𝐒_0`$. And this is the only term that produces the Schiff moment of the nucleus with an odd neutron, like 199Hg. The residual interaction $`F`$ has the form

$$F=C\left(f(r)+f^{}(𝝉_1𝝉_2)+g_s(𝝈_1𝝈_2)+g_s^{}(𝝈_1𝝈_2)(𝝉_1𝝉_2)\right)\delta (𝐫_1𝐫_2),$$
(16)
where $`C=`$300 MeV$`fm^3`$, and $`f(r)=f_{ex}+(f_{in}f_{ex})\frac{\rho (r)}{\rho (0)}`$. Note that the angular dependence of the operators $`𝐒`$ and $`\delta 𝐒`$ is completely different. While $`S_\mu Y_{1\mu }(𝐧)`$, the induced part $`\delta S_\mu `$ is a superposition of spin dependent operators $`\sigma _\mu `$ and $`\{\sigma Y_2(𝐧)\}_{1\mu }`$. For this reason, different parts of the interaction (16) contribute into renormalization of $`𝐒`$ and $`\delta 𝐒`$.

For a spherical nucleus we can separate the angular variables and solve the obtained equations in coordinate space. The equations are

$$S^a(r)=S_0^a(r)+_0^{\mathrm{}}A^{ab}(r,r^{})S^b(r^{})𝑑r^{},$$
(17)
where $`a=p,n`$ and $`S_0^p(r)`$ is the radial part of the Schiff moment operator Eq. (1) multiplied by $`r`$. The particle-hole propagator $`A(r,r^{})`$ was calculated by means of the Green functions of radial Schrödinger equation.

$$A(r,r^{})=\frac{C}{3}Tr_2\{(f(r)+f^{}(𝝉_1𝝉_2))\underset{\kappa jl}{}n_\kappa |jl||Y_1||\kappa |^2rR_\kappa (r)\times $$

$$r^{}R_\kappa (r^{})(G_{jl}(r,r^{}|ϵ_\kappa +\omega )+G_{jl}(r,r^{}|ϵ_\kappa \omega )\}.$$
(18)
Similar equations can be written for $`\delta S(r)`$ \- the local part of the induced moment. They differ from Eq. (18) in type of tensor operators and in residual interaction.

### 3.1 Separation of the spurious component

The integral equation (17) must have zero eigenmode related to the center of mass motion. However, in our case the mean field Eq. (6) is not consistent with the interaction Eq. (16) and, in general, we do not have zero energy for the center of mass motion. The situation can be improved using some freedom in the value of the interaction constant $`f_{in}`$ in Eq. (16). We can fix the value by the condition $`\omega _0=0`$ for the lowest isoscalar dipole mode. Since this procedure is numerical the condition $`\omega _0=0`$ cannot be fulfilled exactly but with finite accuracy. For 199Hg we found that the set $`f_{in}=0.3935`$, $`f_{ex}=2.6`$, and $`f^{}=1.07`$ gives for $`\omega _0`$ the value $`\omega _0=0.1`$ KeV which is really small compared to the energy of dipole transitions.

The finite accuracy in determination of the spurious mode brings another problem for solutions of the Eq. (17). The bare operator $`S_0(r)`$ becomes non orthogonal to the spurious mode transition density. This results in admixture of the spurious component to all solutions of the Eq. (17). The spurious component should be subtracted since it can change the solution considerably. The subtraction can be performed using analytical properties of the solution as a function of $`\omega `$. The solution of the inhomogeneous linear integral equation (17) as a function of $`\omega `$ has the first order poles at $`\omega =\pm E_{ex}`$, where $`E_{ex}`$ are the excitation energies of the RPA modes. Thus, at small $`\omega `$ the solution can be presented as

$$S(r|\omega )=\frac{a(r)}{\omega ^2\omega _0^2}+b(r|\omega =0)+O(\frac{\omega ^2}{E_{ex}^2}).$$
(19)
The first term in Eq. (19) is the spurious component contribution. Let us define the following set of integrals in a complex $`\omega `$-plane along a circle with the radius satisfying the conditions $`\omega _0|\omega |E_{ex}`$.

$$J_n(r)=\omega ^nS(r|\omega )\frac{d\omega }{2\pi ı}.$$
(20)
Using Eq. (19) we can perform the integration analytically As a result we obtain

$$b(r|\omega =0)=J_1(r),$$

$$\omega _0^2=J_3(r)/J_1(r).$$
(21)
The numerical integration in Eq. (20) has been performed using 32-points Gauss formula. The integration radius was $`|\omega |=0.1`$ MeV. The obtained field $`b(r|\omega =0)`$ was practically insensitive to the integration radius. The value of $`\omega _0^2`$ was independent on $`r`$ in first 5 or 6 digits. The loss of accuracy was noticeable only at $`|\omega |\omega _0`$.

Fig.1 shows the result for renormalized $`S(r)=b(r|\omega =0)`$ with subtracted spurious component. It is shown by full line. The dashed line shows the unrenormalized component $`S_0(r)`$. The effects of the core polarization are not large. We found that for 209Bi they change the valence proton contribution by $`15\%`$ which is in fair agreement with the estimates of Ref. . The proton component does not contribute to the Schiff moment of 199Hg since the valence nucleon is a neutron. However, due to neutron-proton residual strong interaction some neutron component is produced by the core polarization. This component is shown in Fig.2.

Qualitatively, the behavior of the neutron component inside a nucleus resembles the behavior of the bare operator in Fig.1. However, the magnitude of the neutron component is smaller and the behavior outside of the nucleus is completely different.

## 4 Results for 199Hg

In Table 1 we show the contributions of the neutron component discussed above to the Schiff moment of 199Hg .

The three columns correspond to three isospin channels. The first row is the contribution from $`\delta R_{dir}`$ produced by $`\delta U_{dir}`$ given by Eq. (7). The second row is the contribution from $`\delta R_{exch}`$ produced by $`\delta U_{exch}`$ given by Eq. (9), and the third row is the sum of them.

In the direct contributions the isospin channel $`T=1`$ dominates, as it was mentioned above. The exchange contributions are small, in general. For the isospin channels $`T=0,2`$ they are comparable to the direct contributions just because the llatter’s are suppressed by the factor $`(NZ)/A`$.

The contributions from $`\delta 𝐒`$ are more significant. In Table 2 we list the contributions from the direct and the exchange parts of the weak interaction $`W`$.

In the first two rows the unrenormalized contributions $`\delta 𝐒_0`$ are shown. Again, the exchange contributions are considerably smaller than the direct ones. The contributions from the channels $`T=1`$, and $`T=2`$ are comparable here although $`T=1`$ contribution is still larger. Next two rows show the renormalized contributions $`\delta 𝐒`$. The effect of the core polarization is significant here. The induced moment $`\delta 𝐒`$ includes spin dependent operators, therefore, the spin-spin part of the residual interaction Eq. (16) is responsible for their renormalization. The spin-spin interaction is repulsive with the constants $`g_s=0.63`$, and $`g_s^{}=1.01`$. The repulsion results in decrease of absolute values of the renormalized contributions. Finally, in the last row we show the contribution of the non-local term $`\delta 𝐒_{NL}`$ Eq. (15). It has the exchange origin as well and its contribution is really insignificant. All the contributions can be summarized as follows

$$S=0.0004gg_00.055gg_1+0.009gg_2[e\mathrm{fm}^3].$$
(22)
The obtained value for the Schiff moment in Eq. (22) cannot be compared directly with the previous calculations with the contact interaction Eq. (5). The reason is in different definition of the dimensionless constants $`gg_i`$ in Eq. (4) and $`\eta _{ab}`$ in Eq. (5). To perform the comparison we redefine the constants $`g_i`$

$$g_i=\frac{Gm_\pi ^2}{\sqrt{2}}\stackrel{~}{g}_i.$$
(23)
With this factor the integration over space of the Yukawa function gives $`1`$, exactly as the integration of a $`\delta (𝐫)`$. Introducing this factor we obtain

$$S=(0.01g\stackrel{~}{g}_00.86g\stackrel{~}{g}_1+0.14g\stackrel{~}{g}_2)\times 10^8[e\mathrm{fm}^3].$$
(24)
This value should be compared with $`S=1.4\times 10^8\eta _{np}`$ from Ref. , and $`S1.6\times 10^8\eta _{np}`$ from Ref. . Remembering that $`\eta _{np}g(\stackrel{~}{g}_0+\stackrel{~}{g}_12\stackrel{~}{g}_2)`$ we conclude that the difference between our result and previous calculations is significant for $`T=0`$, and $`T=2`$ channels. Our values are smaller in absolute value. In order to trace the origin of this difference we repeated our calculations using the contact interaction and omitting the core polarization. The contact interaction was obtained by replacing the Yukawa function in Eq. (4) by the delta-function. The result is

$$S=0.96\times 10^8g(\stackrel{~}{g}_0+\stackrel{~}{g}_12\stackrel{~}{g}_2)[e\mathrm{fm}^3].$$
(25)
If we omit completely the effect of the core polarization in calculations with the finite range interaction Eq. (4) then the only nonzero contribution in the first row of the Table 2 becomes

$$S=0.086g(g_0+g_12g_2)[e\mathrm{fm}^3],$$

that corresponds to

$$S=1.35\times 10^8g(\stackrel{~}{g}_0+\stackrel{~}{g}_12\stackrel{~}{g}_2)[e\mathrm{fm}^3].$$
(26)
Comparing Eq. (25) and Eq. (26) we conclude that the effect of finite weak interaction range is not very significant. The main effect bringing the value of the Schiff moment from that in Eq. (26) to the value in Eq. (24) comes from the core polarization.

The last remark concerns the pairing effects. The nucleus 199Hg has 7 neutron holes in the unfilled shell. One of them is fixed in $`p_{1/2}`$ state and 6 others should be distributed among the other states in the shell. From mass differences we found $`\mathrm{\Delta }_n=0.69`$ MeV. This value of the pairing gap is typical for developed pairing. For dipole transitions the transition energy is large compared to $`\mathrm{\Delta }`$ and the pairing effects are small. They were omitted in Eq. (17). For the induced moment the situation is different. There, the transitions with $`\mathrm{\Delta }J=0,1`$ and $`\mathrm{\Delta }L=0,2`$ are responsible for the core polarization. Such transitions exist inside the last unfilled neutron shell and, due to Pauli blocking, they are sensitive to the details of the shell occupation. For a T-odd operator the effects of pairing in the core polarization can be considered in the way used in Ref. . We found, that the main effect of pairing is in fixing the occupation numbers in the upper unfilled neutron shell. As soon as we keep the occupation numbers fixed the results for $`\delta 𝐒`$ are changed within few percents when we put $`\mathrm{\Delta }_n=0`$.

In summary, we calculated the Schiff moment of 199Hg nucleus using finite range weak interaction and considering the core polarization effects. The effects of the finite interaction range are not very significant for the Mercury nucleus. They do not change the order of magnitude of the Schiff moment calculated with the contact interaction. The effects of the core polarization are two-fold. First, they renormalize the bare operator of the Schiff moment producing a small neutron component due to n-p residual interaction. Second, they produce an induced Schiff moment due to P- and T- violation component in the intermediate single particle states. The induced moment is proportional to both the strong residual and the weak interactions. The effects of the core polarization for the Mercury nucleus are large and they have to be accounted in calculations of P- and T-violating effects in the Mercury nucleus.

## Acknowledgments

The authors appreciate discussions with I.B. Kriplovich during this work.


text:
# Theory for structure and bulk-modulus determination

## I Introduction

Calculations of the bulk structure and the bulk elastic properties play an important role in the physics of condensed matter.alchagirov ; kackell ; persson ; kappaalumina ; pre00 ; wach  Bulk calculations help us understand, characterize, and predict mechanical properties of materials in our surroundings, under extreme conditions, as in geological formations and setting,katahara  and for industrial applications.kappaalumina ; japan  Crystalline materials come in many different structures and, in contrast to isotropic materials, the structure description of crystalline materials may in general need multiple lattice parameters and an atomic basis. In this paper we discuss how to determine the equilibrium structure of a (multiparameter) crystalline material while, at the same time, directly determining the bulk modulus and the bulk modulus pressure derivative. We argue and show that for theoretical structure calculations of multiparameter systems this is simpler and more exact than fitting to (semi-)empirical equations of state (EOS) such as, e.g, the Murnaghan or Birch EOS. In particular, with our direct method there is no need to first determine the hydrostatic path of the system. We further discuss how to include the atomic basis in this process in a natural way.

In crystalline materials described by a single lattice parameter (e.g, monatomic cubic phases) the lattice parameter is a simple function of the unit-cell volume, and the equilibrium volume thus uniquely determines the equilibrium structure, i.e., the value of the lattice parameter. This is not the case when multiple lattice parameters characterize the system and a whole range of lattice-parameter values can form the same unit-cell volume. The equilibrium structure of the material must then be found by fitting and minimizing the free energy within the multidimensional space of lattice parameters.kackell ; wach ; comment  Relevant variables describing the atomic basis (e.g., bond lengths or binding angles) may be included among the parameters, and the full set of lattice parameters and internal (atomic basis) parameters are collected into the vector $`𝒙`$, scaled to dimensionless form. The volume of the unit cell $`V(𝒙)`$ depends in a simple way on the values of the lattice parameters describing the unit cell, but not on the internal atomic configuration. Nevertheless, we here treat the external and internal parameters on an equal footing.

From theory bulk calculations the total energy (per unit cell) $`E(𝒙)`$ is found for a number of structures $`𝒙`$. The elastic response of typical hard crystalline materials corresponds to small deviations $`\delta 𝒙=𝒙𝒙^{(0)}`$ of the structural parameters from the equilibrium structure $`𝒙^{(0)}`$. The observation that the total energy forms a natural potential (hyper-)surface in the parameter space of lattice and internal parameters $`𝒙`$, combined with the accuracy of present-day bulk-calculation methods (such as density functional theory, embedded atom methods, or effective medium theory), then makes it possible to fit the corresponding total-energy variation through the multidimensional fit

$$E(𝒙)=k+\frac{1}{2}M_{ij}\delta x_i\delta x_j+\frac{1}{3!}\gamma _{ijl}\delta x_i\delta x_j\delta x_l+O\left(\delta 𝒙\right)^4,$$
(1)
at controlled accuracy. Here $`k`$, $`M`$, and $`\gamma `$ denote zeroth, second, and third-rank tensors of fitting constants. An additional set of fitting constants are the $`𝒙^{(0)}`$ hidden in $`\delta 𝒙=𝒙𝒙^{(0)}`$ The polynomial fit (1) gives a transparent description of the materials-structure energy variation and directly determines the equilibrium structure $`𝒙^{(0)}`$.

In this paper we exploit and use the structure calculation, i.e., the multidimensional polynomial fit of the total energy (1) for an additional and direct determination of the zero-pressure bulk modulus $`B_0=V(𝒙^{(0)})(p/V)|_{𝒙=𝒙^{(0)}}`$ and its pressure derivative, $`B_0^{}=\{V(p/V)\}/p|_{𝒙=𝒙^{(0)}}`$ at zero temperature.

For a general set of structure parameters, $`𝒙`$, we expand the volume around the equilibrium configuration $`𝒙^{(0)}`$ using the gradient $`𝒈=V(𝒙)|_{𝒙=𝒙^{(0)}}`$ and the Hessian $`H=H\left(V(𝒙)\right)|_{𝒙=𝒙^{(0)}}=\left[\left\{^2V(𝒙)/(x_ix_j)\right\}_{ij}\right]_{𝒙=𝒙^{(0)}}`$ of the volume. We note that derivatives of the volume with respect to the internal parameters vanish, by definition. By providing a systematic treatment of the structural changes induced by the pressure $`p=E/V`$ we extract from the minimum of the zero-temperature enthalpy

$$(𝒙,p)=E(𝒙)+pV(𝒙)$$
(2)
both the bulk modulus

$$B_0=\frac{V(𝒙^{(0)})}{𝒈^TM^1𝒈}$$
(3)
and the bulk-modulus pressure derivative

$$B_0^{}=V(𝒙^{(0)})\frac{3𝒈^TM^1HM^1𝒈\gamma _{ijl}\left(M^1𝒈\right)_i\left(M^1𝒈\right)_j\left(M^1𝒈\right)_l}{\left(𝒈^TM^1𝒈\right)^2}1.$$
(4)
The algorithm outlined above can also be applied to the corresponding direct determinations of general harmonicnye  and anharmonic elastic properties.Forthcoming

Our results both enhance the theory understanding of the crystalline mechanical properties and simplify the desired testing of theory calculations as they combine the formal determination of the crystalline structure \[Eq. (1)\] and of the elastic properties \[Eqs. (3) and (4)\]. For example, from Eqs. (1) and (2), we can directly identify which (internal) structure parameters softens the bulk modulus (3) and we may, in turn, strengthen the materials by suitable chemical or structural modification.

For cases like graphite, where the state-of-the-art DFT based on the generalized gradient approximation (GGA) fails to describe the weak physical interlayer binding, the polynomial fit (1) highlights the intrinsic theory challenges:japan  the total energy has no minimum in the space of lattice parameters, the too soft dependence of the total energy on the graphite-layer separation is quite apparent in the fit. At the same time, for strongly bonded crystalline materials the fit (1) produces usable estimates of the equilibrium crystalline structure and shows the strength of GGA-DFT. The equilibrium volume, the crystallographic parameters, and the bulk modulus, describing the material’s resistance to hydrostatic stress, provide simple experimental tests against which we can compare and calibrate our calculations.

Besides the direct relevance of our results for the description of complex materials our calculations of bulk structure and bulk modulus calculations are also of interest for development of pseudo-potential-based density-functional-theory (DFT) methods and for methods using empirical parameters. There, a first and critical test of the pseudopotential or the empirical parameters is whether the calculations predict a correct materials structure, binding, and elastic properties for the relevant equilibrium configuration. Present DFT scriptsdacapo  can automate some pseudopotential testing for simple materials and symmetries, our formal results generalize such testing of theory accuracy to cases when multiple structural parameters determine the elastic properties.

The outline of this paper is as follows. In Section II we discuss the traditional methods of determining the bulk modulus for single- and multiparameter systems. In Section III we derive our expressions for the bulk modulus and the bulk modulus derivative, Eqs. (3) and (4), for the simple one-parameter problem (e.g., mono-atomic fcc or bcc structures), easily generalized to the $`n`$-parameter problem. In Section IV we proceed to illustrate and test the algorithm on a number of mono- and di-atomic materials based on first-principle DFT calculations and comparison to experiments. Comparisons of $`B_0`$ and $`B_0^{}`$, together with the test of the lattice and structure parameters themselves, represent the typical test of materials-theory accuracy. Section V contains the conclusion.

## II Background

A theory determination of the zero-temperature bulk modulus based on either traditional methodsalchagirov ; murnaghan ; birch  or our formal result (3) is straightforward when one single structural parameter (e.g., the lattice parameter $`a`$) defines the crystalline state. This situation applies for monatomic crystals with simple cubic (sc), face-centered cubic (fcc) and body-centered cubic (bcc) symmetries. Here, the unit-cell volume $`V(a)=qa^3`$ uniquely determines the lattice parameter $`a`$ through a dimensionless number $`q`$ which depends on the crystal symmetry ($`q=1`$, $`q=1/4`$, and $`q=1/2`$ for sc, fcc, and bcc lattices, respectively). All which is required are theory calculations of total energies for a range of $`a`$ values to determine both the equilibrium structure $`a_0`$ and the equilibrium volume $`V^{(0)}`$. The total energy per unit cell, $`E(a)`$ (as in Eq. (1)), can then be expressed as a function of the unit-cell volume, $`E(V)`$.

The general approach is illustrated by the example in Figure 1, which shows the total energy as a function of the lattice parameter $`a`$ for the zinc-blende phase of SiC (3C-SiC), as found from DFT calculations. A parabola fit to the 15 data points closest to the equilibrium value of $`a`$ and a fourth-order polynomial fit to all data points are shown. Fits using the traditional Murnaghan equation of state,murnaghan  integrated to give

$$E_{\text{Murn}}(V)=E_0+\frac{B_0V}{B_0^{}}\left[\frac{(V^{(0)}/V)^{B_0^{}}}{B_0^{}1}+1\right]\frac{V^{(0)}B_0}{B_0^{}1}$$
(5)
or the Birch equation of state,birch  integrated to

$`E_{\text{Bir}}(V)`$ $`=`$ $`E_0+{\displaystyle \frac{9}{8}}B_0V^{(0)}\left[\left(V^{(0)}/V\right)^{2/3}1\right]^2`$ (6)
$`+{\displaystyle \frac{9}{16}}B_0V^{(0)}\left\{B_0^{}4\right\}\left[\left(V^{(0)}/V\right)^{2/3}1\right]^3`$
$`+𝒪\left[\left(V^{(0)}/V\right)^{2/3}1\right]^4,`$
yield, to the eye, curves identical to the fourth-order polynomial fit and are not shown separately. In Murnaghan’s and Birch’s equations (5) and (6) the quantities $`B_0`$, $`B_0^{}`$, and $`V^{(0)}`$, and in some cases also the cohesive energy $`E_0`$, are fitted. Other equations of state traditionally used are mentioned in Refs. alchagirov,  and birch, .

The values of the equilibrium lattice parameter $`a_0`$, and of $`B_0`$ and $`B_0^{}`$ obtained from Eqs. (3) and (4) and from the Murnaghan and Birch fits are included in Table 1. For systems described by one lattice parameter Eqs. (5) and (6) give bulk moduli and bulk-modulus derivatives in close agreement with our present direct approach, Eqs. (3) and (4).

We would like to stress that the moduli $`B_0`$ and $`B_0^{}`$ are formally defined as zero-pressure quantities, and in no way depend on finite-pressure behavior beyond the pressure gradient at $`p=0`$. If we are able to sample our theory system in a sufficiently dense grid around the zero-pressure structure the values of $`V_0`$, $`B_0`$, and $`B_0^{}`$ in (3) and (4) are exact,accuracy  and can be related to the corresponding exact determination of the elastic constants. Fits to empirical EOS may yield results of $`V_0`$, $`B_0`$, and $`B_0^{}`$ that are in good agreement with experimental observations, but they do not necessarily constitute the exact quadratic response.

For materials with multiple structure parameters, the procedure of the traditional approaches further becomes quite awkward as it must be supplemented by a separate discussion of how the experimental conditions define the relevant structural constraint at a given volume, the hydrostatic path $`𝒙=𝒙(V)`$. Moreover, cross-correlations on the $`n`$-dimensional energy surface are ignored in traditional fitting procedures. These procedures are basically a one-dimensional fit in the $`n`$-dimensional space and they are thus more subject to numerical noise in the data points than our approaches based on the multidimensional least-squares polynomial fit (1).kackell

A simple multi-parameter case illustrates this point. Figure 2 both describes the energy surface, Eq. (1), fitted through DFT calculations for Co, and emphasizes the general advantages of our direct approach \[Eqs. (3) and (4)\] compared with the traditional bulk-modulus determinations.alchagirov ; murnaghan ; birch  Materials like Co, which has a non-ideal hexagonal close-packed (hcp) structure, graphite with its layered structure, or the polytypes of SiCkackell  or alumina,kappaalumina  have multiple lattice parameters (plus relevant internal degrees of freedom) which are, of course, no longer uniquely specified by the volume but depend on the general materials conditions. The hydrostatic path defines the system when subject to uniform pressure (as relevant for the bulk-modulus measurements). The hydrostatic path is identified in Fig. 2. The traditional approachesmurnaghan ; birch  proceed by implementing this complex constraint in the equation of state, Eq. (5) or (6), to fit the bulk modulus and its derivative. Instead, we present an explicit determination, Eqs. (3) and (4), based directly on the equation of state (1) expressed as a function of the underlying crystalline structural parameters.

## III Derivation

Our direct bulk-modulus evaluations, Eqs. (3) and (4), are the results of using the pressure (instead of the volume) in a formal identification of the hydrostatic path and then invoking a systematic expansion in small pressure for an explicit specification of $`B_0`$ and $`B_0^{}`$. Today DFT and other materials-theory bulk-calculations are done with high accuracy, and we need only to vary the lattice parameters values slightly around the optimal structure to approximate the total-energy curve by an accurate polynomial fit, Eq. (1). The minimum of the corresponding (zero-temperature) enthalpy (2) can thus be used to directly specify the physically correct structural configuration at any given pressure $`p`$. The set of these optimal structure-parameter values, $`𝒙_{\text{hydro}}(p)`$, trace out the hydrostatic path which, when parameterized by $`p`$, is obtained by simply solving the equation

$$(𝒙,p)=E(𝒙)+pV(𝒙)=0.$$
(7)
We obtain a formal expression for the general (pressure dependent) bulk modulus by taking the derivative along this hydrostatic path

$$B=V(𝒙_{\text{hydro}}(p))\left(\frac{V(𝒙_{\text{hydro}}(p))}{p}\right)^1,$$
(8)
and finally we extract the explicit results for the zero-pressure bulk-modulus $`B_0`$ and for its pressure derivative, $`B_0^{}`$.

We illustrate the general derivation by focusing on a one-dimensional parameter space $`x`$, e.g., for the fcc or bcc one-atomic structure. Then the total energy can be fitted by the polynomial

$$E(x)=k+\frac{1}{2}M(xx^{(0)})^2+\frac{1}{3!}\gamma (xx^{(0)})^3+f(xx^{(0)})$$
(9)
where $`f(xx^{(0)})=𝒪(xx^{(0)})^4`$ contains higher order terms. The coefficients $`k`$, $`M`$, $`\gamma `$, the coefficients of $`f(xx^{(0)})`$, as well as the optimal value of the lattice parameter at zero pressure, $`x^{(0)}`$, are the fitting parameters to be specified, for example, by a set of accurate underlying DFT calculations.

At small pressures, i.e., for lattice and internal parameter values close to the zero-pressure optimal values, the $`pV`$ term in the enthalpy is small and can be regarded as a perturbation of the system. To proceed we introduce a small, non-dimensional and real parameter $`\lambda `$ such that we can write the small pressure as $`p=\lambda p^{(1)}`$ and the lattice-parameter variable as $`x=x^{(0)}+\lambda x^{(1)}+\lambda ^2x^{(2)}+𝒪(\lambda ^3)`$. The variables $`x^{(1)}`$, $`x^{(2)}`$, … are unknown, are functions of the pressure and must be found in the following.

The bulk modulus expression requires calculation of the volume $`V(x)`$ and its pressure derivative. We write the volume in a Taylor-expansion around the zero-pressure solution $`x^{(0)}`$ as

$`V(x)`$ $`=`$ $`V(x^{(0)})+\lambda gx^{(1)}`$ (10)
$`+\lambda ^2\left(x^{(2)}g+{\displaystyle \frac{1}{2!}}\left(x^{(1)}\right)^2H\right)+𝒪(\lambda ^3)`$
with $`g=dV/dx|_{x=x^{(0)}}`$ and $`H=d^2V/dx^2|_{x=x^{(0)}}`$. Here, the pressure dependence enters through the variables $`x^{(1)}`$, $`x^{(2)}`$, …. The pressure derivative of the volume is thus

$$\frac{V(x)}{p}=g\frac{x^{(1)}}{p^{(1)}}+\lambda \left(g\frac{x^{(2)}}{p^{(1)}}+\frac{x^{(1)}}{p^{(1)}}Hx^{(1)}\right)+𝒪(\lambda ^2)$$
(11)
and we determine the variables $`x^{(1)}`$, $`x^{(2)}`$, … by solving the condition on the enthalpy given by (7)

$`0`$ $`=`$ $`\lambda \left(Mx^{(1)}+p^{(1)}g\right)`$
$`+\lambda ^2\left(Mx^{(2)}+{\displaystyle \frac{1}{2}}\gamma \left(x^{(1)}\right)^2+p^{(1)}Hx^{(1)}\right)+𝒪(\lambda ^3).`$
The identity (LABEL:eq:idlambda) must hold for every order and we thus obtain a formal pressure dependence of the lattice parameter

$`x^{(1)}`$ $`=`$ $`p^{(1)}M^1g`$ (13)
$`x^{(2)}`$ $`=`$ $`\left(p^{(1)}\right)^2\left\{M^1HM^1g{\displaystyle \frac{1}{2}}M^1\gamma M^1gM^1g\right\}.`$
Finally, introducing these solutions into (11) we find for $`\lambda =0`$ the isothermal zero-pressure bulk modulus

$$B_0=V(x^{(0)})\left(\frac{V}{p}\right)_{x=x^{(0)}}^1=\frac{V(x^{(0)})}{gM^1g}$$
(15)
and taking the derivative of $`V\left(V/p\right)^1`$ with respect to $`p=\lambda p^{(1)}`$ we find at $`\lambda =0`$

$$B_0^{}=V(x^{(0)})\frac{3gM^1HM^1g\gamma M^1gM^1gM^1g}{\left(gM^1g\right)^2}1$$
(16)
in the case when one (lattice) parameter suffices to describe the unit cell and its atom basis.

The above derivation is straightforwardly generalized to materials systems in which $`n`$ independent lattice and internal parameters determine the structure and the bulk moduli Eqs. (3) and (4), which is our main result. We stress that $`B_0`$ and $`B_0^{}`$ are evaluated at zero pressure and thus the results are exact in spite of the perturbation. $`B_0`$ and $`B_0^{}`$ depend directly on the second order, respectively on the second and third order, coefficients of the energy fit ($`M`$ and $`\gamma `$). We observe that the coefficients of $`f(xx^{(0)})`$ do not enter the expression for the bulk modulus (3) or the pressure derivative of the bulk modulus (4). However, their presence may improve the fit (9), and thereby affect also the coefficients $`M`$, and $`\gamma `$, and thus $`B_0`$ and $`B_0^{}`$. Internal parameters, which describe the positions of the atoms within the unit cell, naturally do not enter the expression of the volume, and thus not the volume derivatives $`𝒈`$ and $`H`$ either, but do affect $`B_0`$ and $`B_0^{}`$ through $`M^1`$.

Higher pressure derivatives of the bulk modulus may be found by taking into account the higher orders of $`\lambda `$ in the Taylor expansions of $`𝒙`$ and the volume. The pressure derivatives will depend on successively higher orders in the polynomial fit. The derivation is straightforward if somewhat tedious.

## IV Examples of applications

As an example of the use of the algorithm for determining $`B_0`$ and $`B_0^{}`$ we evaluate the structure and bulk modulus of a selection of one and two-species materials. We fit data obtained from DFT calculations, described in further detail below, to fourth-order polynomials of the form (1) in $`n`$-dimensional space, where $`n=1`$, 2, 3, or 5.

The pseudopotentials used in DFT calculations may be optimized for various purposes, but should generally yield consistent and transferable accuracy and results. Here, we have used some of the pre-defined pseudopotentialspseudolist  of the open-source DFT program DACAPO.dacapo  The values that we find for the lattice constants, for $`B_0`$, and for $`B_0^{}`$, are collected in Table 1. For reference, the experimental values are also included, as well as the bulk modulii from a Murnaghan and Birch fit along the hydrostatic path.

We have calculated the structure and bulk modulus for three multiparameter systems, as well as for a number of related one-parameter systems: the two-parameter hcp phase of Co, the three-parameter (one internal) wurtzite phase of SiC (2H), the five-parameter (three internal parameters) hexagonal 4H-polytype of SiC, and the one-parameter bcc and fcc phases of Co, the zinc-blende (3C) phase of SiC and the diamond phases of C and Si.

For the DFT calculations we used the plane-wave DACAPO codedacapo  with GGA. For the calculations of the 2H and 4H-polytype of SiC we used $`8\times 8\times 8`$ and $`8\times 8\times 4`$ $`k`$-points, respectively, to describe the Brillouin zone. For all other calculations $`10\times 10\times 10`$ $`k`$-points were used. A uniform energy cut-off of 400 eV, and a conservative choice of fast-Fourier transform (FFT) grid was used. For each evaluation of the optimal structure we calculated by DFT a number of data points for the lattice parameter(s) approximately within $`\pm 10`$% from the expected optimal value(s) of the lattice parameter(s). In the one-parameter systems we calculated 20-30 data points, for the Co hcp-structure 120 data points, for the 2H-polytype of SiC 140 data points, and for 4H-SiC we calculated 7500 data points. The Co calculations were spin-polarized, yielding realistic values of the spin polarization over the range of lattice parameter values considered here.

The possibility of treating the external and internal parameters collectively is important. For example, the variation of internal bond length with pressure might be as important for the total energy as the change in (external) lattice parameters. Often, relaxation of the internal parameters is done with a steepest-descent (or similar) search minimizing the Hellmann-Feynman forces to a certain cut-off at fixed lattice parameters. Although in practice the atomic relaxation will often be a convenient way of obtaining the optimal position of the atoms within the unit cell, the approach has two shortcomings: it introduces a random residual lattice strain, which in turn affects the total energy, and further, the Hellmann-Feynman forces have a non-trivial dependence on the pressure acting on the unit cell and therefore a constant cut-off on the force will not correspond to a constant accuracy of the total energy with varying pressure. Thus a better accuracy — and a consistent choice of accuracy — can be obtained by treating the lattice and internal parameters on an equal footing. This is here done for the 2H and the 4H polytype of SiC. The results are shown in Table 1. For the Murnaghan and Birch values of $`B_0`$ and $`B_0^{}`$ we need to explicitly calculate the hydrostatic path \[in $`(a,c/a,u)`$ and $`(a,c/a,u_1,u_2,u_3)`$ space\] before obtaining the fit. In contrast, we stress that when using (3) and (4) there is no need to explicitly calculate the hydrostatic path. This is here done purely for illustrational purposes.

## V Conclusion

In summary, we have presented a new direct algorithm for a combined determination of structure and bulk moduli $`B_0`$ and $`B_0^{}`$. The lattice constants of a multiparameter system are best found in a least-squares polynomial fit, as previously noticed for the SiC hexagonal polytypes.kackell  We show (a) how to exploit this polynomial fit for a direct determination of the zero-pressure bulk modulus $`B_0`$ and its pressure derivative $`B_0^{}`$, avoiding the calculation of the hydrostatic path and the subsequent one-dimensional fit to this path. We further show (b) how to consistently include internal parameters, such as bond lengths or bonding angles, in the formalism along with the external lattice parameters. In addition, we have evaluated these formal results in explicit cases within our approach, based on DFT calculations.

###### Acknowledgements.

We thank Per Hyldgaard for useful discussions. This project was partly supported by the EU Human potential research training network ATOMCAD under contract number HPRN-CT-1999-00048. Further, E.S. thanks the Swedish Research Council (VR) and the foundations Trygger and Lundgren for economic support.