text:
# Entanglement generated between a single atom and a laser pulse

## I Introduction

The electromagnetic field describing a laser pulse is intrinsically quantum mechanical and may become entangled with an atom when the two systems interact. The consequences of this fact for experiments which employ lasers as control fields has been the topic of several recent papers vanEnk2001 ; Banacloche2002 ; Silberfarb2003 . All agree that these quantum effects are negligible compared to errors induced by spontaneous emission out of the beam. In particular, in our previous work Silberfarb2003 , we calculated the rate at which information is gained about a two-level atom when one measures a near resonance laser pulse that has interacted with it. This measurement strength could be interpreted as the rate at which photons are spontaneously emitted into the modes defined by the laser beam – a small fraction of the total emission rate. Indeed, as commented by Itano Itano2003 , the quantum fluctuations in a laser beam are equivalent to those in the vacuum, and thus using the standard unitary Mollow transformation Mollow1975 , all quantum effects related to the atom-laser interaction in free space are derived from spontaneous emission in the transformed frame.

Though not an important source of error in quantum operations, entanglement is nonetheless generated between the atom and the laser beam. Photons spontaneously emitted into the laser modes change the quantum statistics of those modes in a way that is correlated with the state of the atom. This entanglement might prove useful as a resource in some quantum communication or quantum computation schemeNielsen2000 . Moreover, by determing how much entanglement one can generate in free space, we obtain a benchmark for comparison with other entanglement generation schemes. For example, the entanglement generated in a cavity QED geometry Kimble1998  could be compared to the entanglement that would be generated in the absence of the resonator.

Quantifying the entanglement generated between an atom and a laser in free space is complicated by the multimode nature of the field in a traveling wave geometry. The Hilbert space describing the laser modes is the tensor product of many independent subsystems. Additionally, the laser and atom do not form a closed system, as the other vacuum modes of the field interact with the atom. Since spontaneous emission into these unoccupied modes dominates over quantum features within the laser modes, it is generally inconsistent to neglect their effects when calculating atom-laser entanglement.

To analyze this problem, we use the formalism developed in our previous paper Silberfarb2003 . We review the salient features in the next section. We then proceed to examine the entanglement generated between the laser pulse mode and the atom. This result is compared to similar results obtained using a closed-system single mode picture. We follow this with a calculation of an upper bound on the total entanglement generated between the atom and all modes propagating with the laser beam. A brief summary of results is presented in the concluding section.

## II Paraxial formalism

In many problems it is useful to have a local description of the quantized field. For example, when the field interacts with a well localized atom or when an experiment uses a detector to continuously measure the field, a delocalized plane wave description is awkward. In systems where a physical quantization volume is defined, such as in cavity QED Kimble1998 , a local set of quantized modes is naturally defined. However, in free space no such intrinsic volume is present. A quantization procedure which considers local temporal modes rather than delocalized spatial modes Blow1990  has a limited regime of validity. Instead, through an appropriate coarse-graining of the standard modal decomposition, we can achieve the desired local picture Silberfarb2003 . A brief review of the pertinent points follows.

The procedure begins by separating the field modes into paraxial and nonparaxial modes. In the canonical example a laser beam impinges upon an atom, interacts with it, then illuminates a detector. This optical system defines a set of modes traveling essentially in one direction with a small diffraction angle. We quantize these paraxial modes by associating a Fock space to every member of a coarse-grained set of wave packets of width $`\mathrm{\Delta }z`$ and transverse area $`A`$. These form an approximate complete set for the paraxial subspace if the dynamics on a time scale smaller than $`\mathrm{\Delta }t=c/\mathrm{\Delta }z`$ are negligible. The coarse-graining length is arbitrarily chosen for convenience and generally should not appear in any experimentally relevant result. It may be chosen to be the width of a trapped atom or to correspond to the detector bandwidth $`\mathrm{\Delta }\omega `$, using $`\mathrm{\Delta }z=c/\mathrm{\Delta }\omega `$. This latter definition is natural when considering the entanglement resource than can be detected. We introduce coarse-grained spatially localized mode functions $`\varphi _i(z)`$ and an associated set of annihilation operators $`\widehat{a}_i`$. The Hamiltonian for the paraxial field and atom in the interaction picture and under the rotating wave approximation is,

$`\widehat{H}_{\text{AF}}(t)=\mathrm{}g{\displaystyle \underset{i}{}}\mathrm{\Theta }_i(t)\left(\widehat{a}_i\widehat{\sigma }_++\widehat{a}_i^{}\widehat{\sigma }_{}\right)`$ (1)
$`\text{where }\mathrm{}g=dE_{\text{vac}}=d\sqrt{{\displaystyle \frac{2\pi \mathrm{}\omega _0}{A\mathrm{\Delta }z}}}.`$
Here $`\mathrm{\Theta }_i(t)`$ are window functions selecting out the mode currently interacting with the atom

$$\mathrm{\Theta }_i(t)=\{\begin{array}{cc}1,\hfill & (i1)\mathrm{\Delta }t<ti\mathrm{\Delta }t\hfill \\ 0,\hfill & \text{otherwise}\hfill \end{array}$$
(2)
In this picture, each wave packet mode propagates past the atom, interacts with it according to a Jaynes-Cummings Hamiltonian Jaynes1963 , then flies away. These dynamics are intrinsically irreversible, in contrast to a cavity QED scenario in which excitaitons oscillate between the field and atom reversibly Kimble1998 . The information contained in these modes can act to decohere the atom at a rate $`\kappa \mathrm{\Gamma }\sigma _{\text{eff}}/A`$, where $`\mathrm{\Gamma }`$ is the total spontaneous emission rate and $`\sigma _{\text{eff}}=3\pi /2k_0^2`$ is the effective scattering cross section on resonance. The “measurement strength”, $`\kappa `$, represents the rate at which photons are spontaneously emitted into the paraxial modes.

The nonparaxial modes which constitute the rest of the field are quantized using the standard free space quantization procedure, associating a Fock space to every plane wave. In the following, these modes will be treated as a zero temperature reservoir into which the atom can decay. The only difference from standard spontaneous emission will be that the rate of decay is slightly reduced by the removal of the paraxial modes.

$$\mathrm{\Gamma }^{}=\mathrm{\Gamma }\left(1\frac{\mathrm{\Omega }_p}{4\pi }\right)\mathrm{\Gamma }\kappa ,$$
(3)
with $`\mathrm{\Omega }_p`$ the solid angle subtended by the paraxial beam. Including the nonparaxial modes, the full master equation for the atom-paraxial field system becomes

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{H}_{\text{AF}},\widehat{\rho }]\frac{\mathrm{\Gamma }^{}}{2}\{\widehat{\sigma }_+\widehat{\sigma }_{},\widehat{\rho }\}+\mathrm{\Gamma }^{}\widehat{\sigma }_{}\widehat{\rho }\widehat{\sigma }_+,$$
(4)
where $`\widehat{H}_{\text{AF}}`$ is the Hamiltonian that couples the paraxial field modes to the atom, given in Eq. (1).

In our canonical example the paraxial modes are taken to be in an initial coherent state, corresponding to the quantum state of a laser beam. We perform a Mollow transformation Mollow1975  on the above master equation to remove the initial coherent field amplitude of the laser and transfer it into the Hamiltonian. In the “Mollow picture”, the initial state of the field is the vacuum and the master equation becomes

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{H}_{\text{coh}}+\widehat{H}_{\text{AF}},\widehat{\rho }]\frac{\mathrm{\Gamma }^{}}{2}\{\widehat{\sigma }_+\widehat{\sigma }_{},\widehat{\rho }\}+\mathrm{\Gamma }^{}\widehat{\sigma }_{}\widehat{\rho }\widehat{\sigma }_+,$$
(5)
where $`\widehat{H}_{\text{coh}}`$ is the coherent semi-classical Rabi-flopping Hamiltonian

$$\widehat{H}_{\text{coh}}=\mathrm{}g\alpha \left(\widehat{\sigma }_++\widehat{\sigma }_{}\right).$$
(6)
The Mollow transformation is achieved by a local unitary operator acting on the field alone and so will not affect the atom-field entanglement. Equation (5) thus serves as our of basic dynamical description for entanglement calculations.

## III Quantifying Entanglement

In reviewing the paraxial formalism in the last section we made a distinction between the paraxial and nonparaxial field modes, the paraxial modes being those modes which travel past the atom and impinge upon the detector. We then traced out the nonparaxial modes to recover the master equation (5). The remaining paraxial subsystem is still a very large dimensional Hilbert space, composed of a tensor product of Fock spaces for each coarse-grained wavepacket in the pulse. Such a large Hilbert space has many subsystem decompositions, any one of which could be of interest when calculating entanglement (Fig. 1). In this section we consider two possible decompositions: (i) entanglement between the atom and the mode defined by the laser pulse, (ii) entanglement between the atom and the entire paraxial subsystem. The former is natural to consider when treating the laser pulse as a single mode, analogous to a Jaynes-Cummings cavity QED geometery Kimble1998 . The latter is of interest when relating the atom-laser entanglement to the measurement strength of the laser and to the error in quantum logic induced by a control pulse. In addition, when calculated using this second division, we can find an upper bound on the atom-laser entanglement that can be used as resource for quantum information processing.

To quantify entanglement, we use a monotone known as the “tangle” which can be related to the square of the more familiar concurrence and also to the entanglement of formation Wootters1998 ; Osborne2002 . For overall pure states where at least one of the systems has only two levels (is a qubit) the tangle is equal to the normalized linear entropy (or purity) of the marginal state of either subsystem Rungta2001 ,

$$T\left(|\psi \psi |\right)=2\left(1\text{Tr}_A\left[\text{Tr}_B\left[|\psi \psi |\right]^2\right]\right).$$
(7)
This formula can be generalized to mixed states by taking the convex roof Uhlmann1998

$$T\left(\widehat{\rho }\right)=\underset{|\psi _i,p_i}{\mathrm{min}}\underset{i}{}p_iT\left(|\psi _i\psi _i|\right)\widehat{\rho }=\underset{i}{}p_i|\psi _i\psi _i|.$$
(8)
We choose to use the tangle instead of the full Von Neumman entropy due to its ease of calculation.

## IV Symmetric mode entanglement

### IV.1 Calculation

We consider first the tangle between the atom an the “symmetric” paraxial mode (Fig. 1), defined by the annihilation operator,

$$\widehat{a}_+=\frac{1}{\sqrt{N}}\underset{i=1}{\overset{N}{}}\widehat{a}_i.$$
(9)
This mode defines a laser pulse of duration $`\tau `$. In calculating the tangle in this decomposition all of the other field modes, nonparaxial and nonsymmetric paraxial modes, act as a reservoir which, when traced over, lead to a mixed-state description of our bipartite system. We assume $`\kappa \tau 1`$, implying that the probability for the atom to emit a photon into one of the paraxial modes is very small. In the Mollow picture (used throughout) the Hilbert space of the quantum field may thus be truncated, allowing at most one excitation in the paraxial subsystem. Note that the atom is still allowed to spontaneously emit an arbitrary number of photons into the nonparaxial field modes. Under these assumptions the atom and the symmetric field mode both behave as two-level systems, coupled together in an overall mixed state. The tangle for such two-level systems can be calculated using Wootters’ general formula for two qubits Wootters1998 .

We begin by calculating the atom-symmetric mode density matrix from Eq. (5) for $`\kappa \tau 1`$. In principle, one could solve for the state of the entire atomic-paraxial system, and then trace out the nonsymmetric modes. However, even limiting all paraxial modes to at most one total excitation there are still $`N+1`$ possible field states, $`N`$ being the number of coarse-grained modes. A more sensible approach is to keep only the “dynamic symmetric mode” at any instant (as defined below), tracing out any non-symmetric modes along the way.

Consider first a single time step $`\mathrm{\Delta }t=\mathrm{\Delta }z/c`$, with $`\mathrm{\Delta }z`$ the coarse-graining length. At the start of this interval the system Hilbert space has dimension $`22`$. The system consists of the two-level atom and a quantized symmetric mode of the field composed only of those coarse-grained modes which have already passed the atom, and restricted to one possible excitation. Assuming $`n<N`$ modes have passed, the annihilation operator for this symmetric field mode will be

$$\widehat{a}_{\text{+,n}}=\frac{1}{\sqrt{n}}\underset{i=1}{\overset{n}{}}\widehat{a}_i.$$
(10)
After evolution over a time interval $`\mathrm{\Delta }t`$, we choose a new subsystem division so that we again have a $`22`$ system. The second subsystem will now correspond to the symmetric mode over the $`n+1`$ coarse-grained modes that have now passed the atom. The procedure may then be repeated to build up the full evolution, and the tangle can be calculated at each step.

The procedure for evolving the state from time $`t`$ to $`t+\mathrm{\Delta }t`$ is as follows. Suppose that $`n<N`$ time intervals have passed. The density operator $`\widehat{\rho }_n`$ will be spanned by the basis states

$$|e|0_{+,n}|g|0_{+,n}|e|1_{+,n}|g|1_{+,n}.$$
(11)
During the next interval a new coarse-grained mode impinges upon the atom. This mode starts in the vacuum (in the Mollow picture) so the new state for this composite system is

$$\widehat{\rho }=\widehat{\rho }_n|0_{n+1}0|.$$
(12)
Through the next time interval this state evolves under a combination of the Jaynes-Cummings Hamiltonian which couples the atom to the quantized mode, a semi-classical coherent atom-laser interaction, and spontaneous emission into the nonparaxial modes. We model this by a short time integration of the master equation, Eq. (5). During this time the density matrix must be expanded to include the newly interacting mode. The basis states include the previous set, Eq. (11), tensored with $`|0_{n+1}`$ for the newly added mode, plus the two new states

$$|e|0_{+,n}|1_{n+1}|g|0_{+,n}|1_{n+1}.$$
(13)
States of the form $`|1_{+,n}|1_{n+1}`$ are higher order terms that are ignored under the assumption that $`\kappa \tau 1`$. This gives a total of six basis states for the density operator. After performing the evolution over this time interval we trace out the nonsymmetric modes of the extended paraxial subsytem. To do this, we perform a unitary transformation on the field operators to create a new tensor product decomposition. Instead of the decomposition described in Eqs. (11,13) consisting of the $`n^{\text{th}}`$ symmetric mode $``$ the $`(n+1)^{\text{th}}`$ coarse-grained mode, we choose symmetric and antisymmetric modes described by the annihilation operators

$$\widehat{a}_{\pm ,n+1}=\sqrt{\frac{n}{n+1}}\widehat{a}_{+,n}\pm \frac{1}{\sqrt{n}}\widehat{a}_{n+1}.$$
(14)
The transformation which achieves this is

$$U=\left(\begin{array}{ccc}1& 0& 0\\ 0& \sqrt{\frac{n}{n+1}}& \sqrt{\frac{1}{n+1}}\\ 0& \sqrt{\frac{1}{n+1}}& \sqrt{\frac{n}{n+1}}\end{array}\right).$$
(15)
Finally $`\widehat{\rho }_{n+1}`$ can be obtained by tracing out the antisymmetric field mode.

The mixed state entanglement between the atom and the symmetric field mode can be calculated using Wootters formula for the tangle Wootters1998  directly from this density matrix. The results will depend on several physical parameters: the area of the beam, the intensity of the laser $`|\alpha |^2`$, and the duration of the pulse. Additionally, the numerical accuracy will depend upon the length of the time step used in the calculations, or equivalently the number of time steps $`N`$. To maintain the paraxial approximation, we choose the area of the beam to be $`\overline{A}=A/\sigma _{\text{eff}}=1000`$. This corresponds roughly to the quadrupole transition $`S_{1/2}`$ to $`D_{5/2}`$ in $`{}_{}{}^{40}\text{Ca}_{}^{+}`$ considered by van Enk and Kimble vanEnk2001 , for which a laser focused to a spot size of 100 $`\mu \text{m}^2`$ corresponds to $`\overline{A}1600`$. For simplicity we consider only $`\pi /2`$ pulses, such that $`g\alpha \tau =\pi /2`$. With the relations $`\kappa =\mathrm{\Gamma }/\overline{A}`$ and $`g=\sqrt{\kappa /\mathrm{\Delta }t}=\sqrt{\kappa N/\tau }`$, assuming a fixed $`N`$, we have a direct relationship between $`\alpha `$ and $`\tau `$. We choose as our free parameter the emission probability, $`\mathrm{\Gamma }\tau `$, which sets the degree of entanglement between the atom and laser.

A plot of the tangle versus $`\mathrm{\Gamma }\tau `$ is given in Fig. 2. Each point in the plot corresponds to the total tangle measured just after the $`\pi /2`$ pulse has completely passed the atom. For very short pulses, the growth in the tangle is linear with $`\mathrm{\Gamma }\tau `$ for all initial conditions, leading to a unit slope on the log-log plot. However the constant of proportionality is different. This difference can be attributed to the probability of the atom being in the excited state during the pulse interval – the more excitation the higher the rate of growth. With a greater probability of spontaneous emission, the behavior deviates strongly from this linear trend, peaking roughly where the pulse length and decoherence time are comparable. After this time the entanglement drops precipitously, and for at least one initial condition, the atom and symmetric field mode become separable. The maximum emission parameter shown in the plot is $`\mathrm{\Gamma }\tau =\sqrt{\overline{A}}=\sqrt{1000}`$, which corresponds to $`\kappa \tau =1/\sqrt{\overline{A}}=1/\sqrt{1000}`$; for a larger emission probability, our truncation of the field Hilbert space to one excitation breaks down.

### IV.2 Approximations

The calculation in Sec. IV.1 yields the entanglement between the atom and the laser mode as it propagates past the atom. It is a complex procedure because we have explicitly separated the two sources of decoherence: spontaneous emission into the nonparaxial modes and spontaneous emission into the paraxial but nonsymmetric modes. We did this because spontaneous emission into the symmetric-paraxial mode is the source of the entanglement and we must be careful not neglect this important effect. Nonetheless, emission into the nonsymmetric-paraxial modes is a small fraction of the total loss since $`\kappa \mathrm{\Gamma }^{}`$. We can consistently maintain the coherent coupling between the atom and the symmetric mode while neglecting emission into the nonsymmetric mode since it acts only as a small source a decoherence, a process dominated by off-axis emission. We thus lump all sources of decoherence together, with total rate $`\mathrm{\Gamma }=\kappa +\mathrm{\Gamma }^{}`$. In this case the Master equation becomes,

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{H}^{},\widehat{\rho }]+\frac{\mathrm{\Gamma }}{2}\left(2\widehat{\sigma }_{}\widehat{\rho }\widehat{\sigma }_+\widehat{\sigma }_+\widehat{\sigma }_{}\widehat{\rho }\widehat{\rho }\widehat{\sigma }_{}\widehat{\sigma }_+\right),$$
(16)
where $`\widehat{H}^{}=\mathrm{}\mathrm{\Omega }(\widehat{\sigma }_++\widehat{\sigma }_{})/2\mathrm{}g(\widehat{a}\widehat{\sigma }_++\widehat{a}^{}\widehat{\sigma }_{})`$ is the driven Jaynes-Cummings Hamiltonian. This equation describes a cavity QED dynamic plus decay into a reservoir Kimble1998 . When the evolution of the tangle is calculated using this equation of motion, the results are in complete agreement with those obtained in Sec. IV.1 to within numerical accuracy.

A more radical simplification is to ignore the decoherence term in Eq. (16) altogether when the interaction time is small compared to the atomic lifetime. The remaining dynamics represents a closed system, with the atom and the symmetric mode evolving unitarily.

$$\frac{\widehat{\rho }}{t}=\frac{i}{\mathrm{}}[\widehat{H}^{},\widehat{\rho }]$$
(17)
For this evolution, the initial bipartite pure state remains pure and the entanglement can be calculated in a straightforward manner. Such an approximation is routinely made for atoms in cavities in the strong coupling regime, where $`g\mathrm{\Gamma }`$ Kimble1998 . In free space we can express $`g`$ in terms of $`\kappa `$; the strong coupling approximation requires $`\kappa \mathrm{\Gamma }^2\tau `$, or equivalently $`\mathrm{\Gamma }\tau \overline{A}1`$. The beam area may be chosen arbitrarily large, however, so that this inequality need not hold. One is not in the strong coupling regime in free space and so agreement between the entanglement calculated using closed and open system models should not be expected apriori. This has been the source of some recent controversy Itano2003 ; vanEnk2003 ; Banacloche2003 .

To compare the numerical results obtained using Eq. (5) with those obtained using the closed system evolution Eq. (17) we first calculate the tangle analytically using the closed system dynamics. Though one could measure the closed system entanglement using the von Neumann entropy of the marginal state as was done in vanEnk2001 , for direct comparison with our results, we use the tangle given in Eq. (7). A simple calculation yields,

$$T=\kappa \tau \frac{\left(2\sigma _y+\pi \right)^2+\sigma _x^2\left(4\pi ^2\right)}{\pi ^2}+O\left((\kappa \tau )^2\right),$$
(18)
where the Pauli matricies describe the two-level atom in the usual way. Expectation values in this equation are to be taken using the initial state of the atom. The results of this approximation are shown in Fig. 2. For the parameters where all approximations are valid, we see a simple linear growth of the entanglement with emission probability. Though slightly larger, the entanglement calculated via the closed system dynamics is in very good agreement, $`O(\mathrm{\Gamma }\tau )`$, with our numerical results for the full atom-symmetric field mode tangle when $`\mathrm{\Gamma }\tau 1`$, as seen in Fig. 2. This agreement holds even when the system is not in the strong coupling regime, $`\overline{A}\mathrm{\Gamma }\tau >1`$.

To understand why the entanglement calculated using closed and open system dynamics agree so well even when the strong coupling assumption is violated, consider the processes which lead to the loss of entanglement. The inherent dynamics of the atom-laser mode system can be represented by a bipartite system, with subsystems $`A,B`$, coupled to a reservoir of other electromagnetic modes, $`R`$, as depicted in Fig. 3. System $`A`$ is coupled to system $`B`$ with strength $`\kappa `$, while it is coupled to the reservoir $`R`$ with strength $`\mathrm{\Gamma }`$. System $`B`$ is not coupled directly to the reservoir, coupling only indirectly through system $`A`$. The total system consisting of both subsystems $`A,B`$ will start in some pure state. From the perspective of quantum trajectories Carmichael1993 , there are two possible ways in which the reservoir $`R`$ can reduce the entanglement generated by the dynamics between $`A`$ and $`B`$. In the first method the nonunitary but deterministic no-jump evolution of system $`A`$ leads to a reduction of this entanglement with probability $`1e^{\mathrm{\Gamma }\tau }`$. The total loss of entanglement is thus proportional to $`\mathrm{\Gamma }\tau T`$ for $`\mathrm{\Gamma }\tau 1`$, where $`T`$ is the tangle generated in the closed system dynamics between $`A`$ and $`B`$. Since entanglement will be generated in the closed system at a rate proportional to the coupling $`\kappa `$, the total loss of entanglement is of order $`O(\mathrm{\Gamma }\tau \kappa \tau )`$. Alternatively, a nondeterministic jump in a quantum trajectory can reduce entanglement if it leads to a particular statistical mixture of states as described below. The probability of this particular type of correlated quantum jump is, however, limited by the smaller decoherence rate of the two subsystems $`A,B`$. In this case system $`B`$ can only decay through system $`A`$ in the interval $`\tau `$ at a rate proportional to $`\mathrm{\Gamma }\tau \kappa `$. The amount of entanglement that can be destroyed through this process will then be of order $`O(\mathrm{\Gamma }\tau \kappa \tau )`$. Thus both deterministic decay and quantum jumps lead to entanglement loss of order $`\mathrm{\Gamma }\tau \kappa \tau `$. This is sufficiently small to ensure that the closed system dynamics generates the same entanglement as the open system dynamics when $`\mathrm{\Gamma }\tau 1`$ and $`\kappa \tau 1`$, independent of the relation between them, i.e. even when the strong coupling approximation breaks down.

The above explanation can be made rigorous for the atom-symmetric field mode system. Assume that $`\mathrm{\Gamma }\tau 1`$ and $`\kappa \tau 1`$, though $`\overline{A}\mathrm{\Gamma }\tau `$ may be large. For $`\mathrm{\Gamma }\tau 1`$ the state of the system after the pulse has passed may be represented by a density operator of the form

$$\widehat{\rho }=(1ϵ)|\psi \psi |+ϵ\widehat{\rho }_{\text{jump}}.$$
(19)
Here $`|\psi `$ is the state evolved under the closed system unitary dynamics, while $`\widehat{\rho }_{\text{jump}}`$ represents the state conditioned on a jump having occurred during the interaction time. The probability for a jump is denoted $`ϵ`$ and will be proportional to $`\mathrm{\Gamma }\tau `$. These assumptions alone are insufficient to ensure that entanglement is preserved. As an example, consider the entangled state $`|\psi =(1\xi /2)|g1+\sqrt{\xi }|e0`$. This state has tangle $`T4\xi `$ for $`\xi 1`$. Note that the state $`|\psi `$ is equivalent under local unitaries to any other state with the same entanglement. If we choose $`\widehat{\rho }_{\text{jump}}=1/2|e1e1\left|+1/2\right|g0g0|`$, the statistical mixture $`\widehat{\rho }`$ given in Eq. (19) is separable for $`\sqrt{\xi }ϵ`$, as can be seen by its positive partial transpose Peres1996 ; Horodecki1996 . For our particular case, an atom in free space, $`\sqrt{\xi }g\tau `$, since entanglement is generated by the atom-field coupling. So for this choice of $`\widehat{\rho }_{\text{jump}}`$, the system can be made separable when $`g<\mathrm{\Gamma }`$, or equivalently $`\overline{A}\mathrm{\Gamma }\tau >1`$.

Since our calculations shows that the open system maintains the same entanglement as the closed system, even when $`\overline{A}\mathrm{\Gamma }\tau >1`$, the dynamics must be restricted such that these special statistical mixtures are rare. To see this, note that the post-jump state we considered above must be generated by a highly correlated jump process. The probability of such a process is very low. Indeed, to good approximation, we expect $`\widehat{\rho }_{\text{jump}}`$ to be of the form

$$\widehat{\rho }_{\text{jump}}=\widehat{\rho }_{\text{atom}}|\text{vac}\text{vac}|$$
(20)
where $`\widehat{\rho }_{\text{atom}}`$ is some arbitrary (possibly mixed) state of the atomic subsystem. This ignores the process of spontaneous emission into a nonparaxial mode followed by subsequent excitation of a paraxial mode. This higher order process may be safely ignored when we are considering entanglement of order $`\kappa \tau `$, and $`\sqrt{\kappa \tau }\kappa \tau \mathrm{\Gamma }\tau `$, which is satisfied when $`\mathrm{\Gamma }\tau 1`$. Inserting Eq. (20) into Eq. (19) the tangle may then be calculated using Wootters’ formula, assuming that $`|\psi `$ is an arbitrary state that can be generated using the closed system dynamic. To lowest order, $`O(\kappa \tau \mathrm{\Gamma }\tau )`$, this agrees with the tangle obtained using closed system dynamics.

### IV.3 Bounding the Total Entanglement

We now consider the entanglement that is generated between the atom and all paraxial modes. This is the physically relevant quantity that determines the error rate in atom-laser control, modulo spontaneous emission into the nonparaxial modes. In addition, it gives a measure of the dynamically generated entanglement resource that is accessible in principle. A calculation of this total entanglement is significantly more complex than the calculation performed in Sec. IV as the Hilbert space associated with all coarse-grained paraxial modes is much larger than that of the single symmetric mode. Even if we assume $`\kappa \tau 1`$, and thereby restrict paraxial field modes to a single excitation as in Sec. IV, the remaining Hilbert space has dimension $`2(N+1)`$. A calculable measure of entanglement for a mixed state of a $`2N`$ system for $`N>2`$ is not known except for some special cases which do not apply here Rungta2003 ; Vollbrecht2001 . We thus resort to calculating a bound for the total entanglement, using the tangle as our measure.

The general difficulty with calculating the entanglement of a mixed state arises because the ensemble decomposition of a density operator into a statistical mixtures of pure states $`\widehat{\rho }=_ip_i|\psi _i\psi _i|`$ is not unique. Averages of a pure state entanglement monotone for different ensemble decompositions do not usually lead to the same result and one must thus perform a difficult minimization search to find the actual entanglement. Since the tangle is a convex function, however, any ensemble decomposition provides an upper bound according to,

$$T(\widehat{\rho })\underset{i}{}p_iT\left(|\psi _i\psi _i|\right)\{\psi _i,p_i\}$$
(21)
Given a known ensemble decomposition, we can calculate the right hand side by averaging the tangle for the given bipartite pure states, over the appropriate probability distribution, using Eq. (7).

A natural ensemble decomposition is obtained by integrating the master equation of the atom-paraxial system using the quantum trajectory method Carmichael1993 . A stochastic Schrödinger equation describes the evolution of an open quantum system with nonunitary evolution due to its coupling to the environment, and conditioned on a fictitious measurement performed on the environment with unit efficiency. For each measurement record $`M(t)`$ there is a corresponding pure state trajectory $`|\psi _M(t)`$. Averaging over the all possible measurement records results in a mixed state description, with density operator given by the convex combination

$$\widehat{\rho }(t)=\underset{M}{}p_M(t)|\psi _M(t)\psi _M(t)|.$$
(22)
From this construction one may immediately write down an upper bound on the tangle using Eq. (21).

Brute force calculation of all the required trajectories is straight forward but prohibitive. If we choose the fictitious measurement to be photodetection, a measurement record will consist a sequence of bits denoting count/no-count decisions in each interval $`(t,t+\mathrm{\Delta }t)`$ (ignoring the possibility of two or more photodetections as it is of order $`\delta t^2`$). Assuming there are $`N=1000`$ such intervals then there are $`2^{1000}`$ trajectories. Our calculation is made tractable, however, by noting that only the most recent jump in any trajectory will contribute to the tangle; the state of any field modes which passed the atom prior to the last jump will not be entangled with the atom since the post-jump state is separable (the atom in the ground state and the field in a coherent state). The evolution after a jump has occurred will be independent of when the jump occurred as the system is Markov, and all of the quantized paraxial modes prior to interaction with the atom are in the same coherent state. It is thus sufficient to calculate two no-jump trajectories – the initial trajectory $`|\psi _A(i)`$ and the post-jump trajectory $`|\psi _B(i)`$. Both represent deterministic evolution according to a nonHermitian Hamiltonian, including the coupling of the system to the nonparaxial modes assuming no quantum jumps into the nonparaxial modes. They differ only in their initial conditions: $`|\psi _A(i)`$ is the evolution given the initial state of the atom at $`t=0`$ and $`|\psi _B(i)`$ is the evolution starting with the atom in the ground state right after a jump. An upper bound on the tangle then follows from Eq. (21),

$$T(\widehat{\rho })\underset{1=1}{\overset{N}{}}P_{LJ}(i)T_{PJ}(Ni)+P_{NJ}T_{NJ},$$
(23)
where $`P_{LJ}(i)`$ and $`P_{NJ}`$ are respectively the probabilitiy that the last jump occurred during the time interval $`i`$ and the probability that no jump occurred during the entire evolution. The tangle associated with no jumps during any interval follows from Eq. (7), $`T_{NJ}=T\left(|\psi _A(N)\right)`$, whereas the post-jump tangle is given by $`T_{PJ}(Ni)=T\left(|\psi _B(Ni)\right)`$.

The relevant probabilities that appear in Eq. (23) may also be calculated using the fiducial trajectories. In particular under nonunitary evolution without normalization, the probability that no jump occurs in a given trajectory up to time $`t`$ is $`\stackrel{~}{\psi }(t)^2`$. Thus, the probability that no jump occurs during any of the $`N`$ intervals is $`P_{NJ}=\stackrel{~}{\psi }_A(N)^2`$. The probability that the last jump occurred in interval $`i`$ is given by,

$$P_{\text{LJ}}(i)=P_{NJ}(N|i)P_J(i),$$
(24)
where $`P_{NJ}(N|i)=\stackrel{~}{\psi }_B(Ni)^2`$ is the conditional probability that no jumps occur between the interval $`i`$ and $`N`$ given that a jump occurred in interval $`i`$. The total probability that a jump occurred at time step $`i`$, $`P_J(i)`$, will satisfy the following equations:

$`P_J(i)`$ $`=P_A(i)+{\displaystyle \underset{j}{}}P_J(i|j)P_A(j),`$ (25a)
$`P_J(i|j)`$ $`=P_B(ij)+{\displaystyle \underset{k>j}{}}P_J(i|k)P_B(kj).`$ (25b)
The first equation expresses the total probability for a jump in interval $`i`$ as the sum of the probability that the system had its first jump in interval $`i`$, $`P_A(i)`$, plus the probability that it first jumped in interval $`j<i`$ then some time later jumped in interval $`i`$. $`P_J(i|j)`$ is thus the conditional probability for a jump at $`i`$ given a previous jump occurred in interval $`j`$. Equation (25b) provides a recursive relation for $`P_J(i|j)`$ in terms of, $`P_B(ij)`$, the probability for a jump to occur $`ij`$ intervals after a previous jump occurred, with no intervening jumps. The probability $`P_A(i)`$ can be easily calculated from the fiducial trajectory $`|\psi _A(i)`$ and similarly $`P_B(ij)`$ can be calculated from $`|\psi _B(ij)`$ Carmichael1993 . Since $`P_J(i|j)`$ is only a function of $`ij`$, as follows from the Markov property, the recursion relation (25b) can be solved independently. We find a suitable solution by truncating the recursion at $`n`$ steps corresponding to at most $`n+1`$ jumps into the nonparaxial modes. Solving for $`P_{LJ}(i)`$, using this truncated solution, we can then obtain an upper bound on the tangle according to Eq. (23).

Our upper bound on the tangle is plotted in Fig. 4 as a function of time. For short times, linear behavior is exhibited, as in the previous calculations for the symmetric tangle. Note, however, that the rate of entanglement generation is always greater than the symmetric calculation. For short pulse durations $`\mathrm{\Gamma }\tau 10^2`$ there is only a $`5\%`$ difference in rates for the initial condition $`|g+i|e`$, the state that is most excited during the interval. All of the plotted upper bounds turn over at roughly the same point, slightly past where the pulse duration equals the decoherence time, and then converge on a single limit for large decoherence probability. This behavior is an artifact our calculation – for long times all information about the initial state has been lost since all trajectories with appreciable probability have undergone multiple jumps.

## V Summary

Entanglement is generated between a laser beam and an atom in free space due to the atom’s spontaneous emission of photons into the paraxial modes, superimposed on the laser’s coherent field amplitude. Through the use of judicius approximations, we were able to quantify this entanglement, a fundamentally hard problem for a large dimensional open quantum system. In particular we calculated the entanglement between the atom and the field mode defined by the laser pulse. Of particular interest is how our measure of a entanglement is reproduced under further simplifying assumptions. We considered two models: a lumped decoherence model and a closed-system model. In the former all sources of decoherence, paraxial and nonparaxial, are lumped into a single decoherence term. This is in excellent agreement with the full calculation when the paraxial approximation holds, $`\overline{A}1`$. This is as expected, since the vast majority of the decoherence is due to the nonparaxial modes. More surprisingly, a closed system model which treats the quantized atom laser-pulse dynamics via a single mode Jaynes-Cummings Hamiltonian is in excellent agreement when $`\mathrm{\Gamma }\tau 1`$, even when the strong coupling assumption was violated. This result was unexpected since the naive picture would argue that, in the absence of the strong atom-mode coupling, $`\overline{A}\mathrm{\Gamma }\tau >1`$, spontaneous emission into other modes should destroy entanglement. We have shown, however, that the particular form of the coupling between the vacuum modes and the atom-laser system in fact preserves the entanglement. Thus, there is no contradiction between the observation that all quantum effects in the atom-laser interaction are due to spontaneous emission Itano2003  and the fact that atom-laser entanglement can be accurately modeled using a single mode Jaynes-Cummings dynamic vanEnk2001 ; Banacloche2002  when $`\mathrm{\Gamma }\tau 1`$. Such agreement is ensured by the fact that spontaneous emission is the only decoherence mechanism.

Finally, we placed an upper bound on the total entanglement that is generated between the atom and the full set of paraxial field modes. We used the quantum trajectory method to calculate a dynamical unravelling of the master equation, which provides a natural ensemble decomposition of the bipartite density matrix. Averaging the pure state tangle over this ensemble leads to an easily calculated upper bound. It remains to be seen how tight this bound is. Nonetheless, it provides a useful benchmark on the total entanglement resource that is available for quantum information processing applications.

###### Acknowledgements.

The authors would like to thank Tracey Tessier and Aldo Delgado-Hilgado for useful discussions and suggestions. This work was supported by NSF Grant No. PHY-009569 and the National Security Agency (NSA) and the Advanced Research and Development Activity (ARDA) under Army Research Office (ARO) Contract No. DAAD19-01-1-0648. AS was supported in part as a QuaCGR Fellow under ARO Contract No. DAAD119-02-1-0213.


format: text

text:
# Bernoulli potential in type-I and weak type-II superconductors: Surface dipole

## I Introduction

The Hall voltage is commonly used to measure the concentration of charge carriers in conductors. In superconductors this method is not applied, because of a missing theoretical support. In this paper we show that the Bernoulli potential, which is closely related to the Hall voltage, can be used to the same end in superconductors.

Any ideal collision-less electric fluid should exhibit a finite Hall voltage, but superconductors seem to escape this theoretical conclusion. The zero Hall voltage was first reported by Kamerlingh Onnes and Hof already in 1914. Though later analyses showed that their samples were in the mixed state what obscures an interpretation of their experiment, the zero Hall voltage was confirmed anyway.

From a theoretical point of view, it was clear that there has to be a voltage balancing the magnetic pressure which acts on electrons via the Lorentz force. The report of Kamerlingh Onnes and Hof thus stimulated various speculations about the missing Hall voltage, see the critical review by Lewis. From various concepts proposed we mention the so-called contact potential – a potential step at the interface – which was expected to cancel the Hall voltage. Such a potential step might exist only if a charge dipole is formed at the interface of a superconductor. The surface dipole we discuss here is a similar concept.

The explanation of the zero Hall voltage turned out to be very simple. By contacts one monitors the differences in the Gibbs chemical potential, often called the electro-chemical potential. The Gibbs potential is composed of three components: the electrostatic potential, the kinetic energy and the correlation energy. None of these components is constant in the presence of diamagnetic currents, but their sum is constant in equilibrium in agreement with the observed zero Hall voltage.

To eliminate the kinetic and correlation energies, Hunt proposed to access the electrostatic potential with a contact-less method called the Kelvin capacitive pickup. The first capacitive measurements appeared soon and they successfully proved the existence of a non-zero electrostatic potential.

Both experiments were done at temperatures well below $`T_c`$, where the electrostatic potential has a simple form resembling the Bernoulli law,

$`e\phi ={\displaystyle \frac{1}{2}}mv^2.`$ (1)
To avoid confusion we note that the electrostatic potential in equilibrium superconductors is called the Bernoulli potential for brevity, even if its actual form does not coincide with the Bernoulli law.

None of the early experimental data were sufficiently accurate to allow for a discussion of possible corrections to the plain Bernoulli potential (1). Nevertheless, the authors made some conclusions in this direction and we find it necessary to comment on them in more detail.

### A Bok and Klein

Bok and Klein claimed that their data agree with the plain Bernoulli potential (1). This conclusion has to be taken with reservations, however, because they measured the electrostatic potential as a function of the magnetic field $`B`$ at the surface. They evaluated the velocity $`𝐯`$ of the superconducting electrons (briefly called the condensate velocity) from the London condition $`m𝐯=e𝐀`$ and the exponential decay $`𝐀=𝐀_0\mathrm{e}^{x/\lambda _0}`$, using $`𝐁=\times 𝐀`$. At low temperatures, the London penetration depth depends on the density $`n`$ of pairable electrons, $`\lambda _0^2=m/(e^2\mu _0n)`$, therefore their experimental result can be expressed in terms of the magnetic pressure, $`en\phi =B^2/(2\mu _0)`$.

According to the above arguments, the relation of the electrostatic potential to the magnetic pressure seems to be a consequence of the plain Bernoulli potential and the London theory. As we show below, the relation to the magnetic pressure is very general and holds also under conditions when neither the plain Bernoulli potential nor the London theory applies.

In fact, Bok and Klein measured on indium which, being a type-I superconductor, is not fully covered by the London theory. Moreover, they swept the magnetic field from zero to the critical value, while the plain Bernoulli potential and the simple form of the London theory used above are restricted to low magnetic fields. At high fields, the condensate density at the surface is suppressed what results in $`v^4`$ and higher-order contributions to the electrostatic potential.

Briefly, Bok and Klein have observed the magnetic pressure. But the link of their experiment to the Bernoulli potential (1) has to be taken with caution.

### B Brown and Morris

Brown and Morris have used a different setup which allowed them to achieve a much higher precision. They did not control the magnetic field but the current in a thin wire. This current was scaled with the critical current. They announced in 1968 that their data reveal about 20% deviations from the screened Bernoulli potential,

$`e\phi ={\displaystyle \frac{n_s}{n}}{\displaystyle \frac{1}{2}}mv^2,`$ (2)
discussed in more detail below.

It should be noted that Brown and Morris expected deviations which were predicted from the BCS theory in the same year. Adkins and Waldram had studied the electrostatic potential from changes of the BCS gap due to a current and they recovered the plain Bernoulli potential for zero temperature, while for finite temperatures they indicated a presence of additional contributions. They were not, however, capable to derive these contributions in an explicit form or to estimate their amplitudes.

Some corrections to the Bernoulli potential (1) were derived already before the BCS studies. Historically the first is the theory of Sorokin from 1949 which covers the majority of effects recovered later. Although this paper is mentioned by London, later it became forgotten. In 1964 van Vijfeijken and Staas took into account that the electrostatic field acts on normal electrons and arrived at the so-called quasiparticle screening.

The quasiparticle screening is represented by the fraction of superconducting electrons $`n_s/n`$ by which eq. (2) differs from the simple Bernoulli law (1). In the experiment of Brown and Morris it accounts for 6% of the observed potential. In spite of its small magnitude, the quasiparticle screening is important with respect to the concept of the magnetic pressure. At finite temperatures one has to take into account that the London penetration depth also depends on the density of condensate, $`\lambda ^2=m/(e^2\mu _0n_s)=\lambda _0^2n/n_s`$. Combining the screened Bernoulli potential (2) with the condensate velocity $`v`$ found from the London theory, one finds that the electrostatic potential is temperature independent, i.e., it is given by the magnetic pressure with the density $`n`$ of pairable electrons. Despite this importance of the quasiparticle screening, it should be noted, however, that the quasiparticle screening is not responsible for the 20% deviations announced by Brown and Morris and discussed above.

Corrections capable to explain the observed potential were first obtained by Rickayzen, who showed that the electrostatic potential includes a pairing contribution,

$`e\phi ={\displaystyle \frac{n_s}{n}}{\displaystyle \frac{1}{2}}mv^2={\displaystyle \frac{n_s}{n}}{\displaystyle \frac{1}{2}}mv^24{\displaystyle \frac{n_n}{n}}{\displaystyle \frac{\mathrm{ln}T_c}{\mathrm{ln}n}}{\displaystyle \frac{1}{2}}mv^2.`$ (3)
The pairing term dominates close to the critical temperature $`T_c`$, because $`n_s=n(1t^4)`$, with $`t=T/T_c`$, while $`n_n=nn_s=nt^4`$.

According to (3), from $`\phi `$ close to $`T_c`$ one may deduce the density dependence of $`T_c`$, which would be very attractive with respect to designing new materials. Indeed, this important material property is otherwise deducible only from measurements applying a hydrostatic pressure or adding impurities to crystals. Unlike the later methods, a measurement of the Hall voltage or, respectively, the Bernoulli potential does not affect the electronic bands, the phonon spectrum or the electron-phonon interaction, therefore it offers a uniquely clear information about the material.

All expectations were chilled by the next paper of Morris and Brown. They admitted that deviations announced in the first paper were due an incorrect estimate of the critical current and presented new highly accurate data for a wide range of temperatures. The observed electrostatic potential is perfectly equal to the magnetic pressure and exhibits no pairing contribution. They reported that this behavior is common to both type-I and weak type-II superconductors and for the magnetic field up to the critical value.

### C Surface dipole

The disagreement between theory and experiment remained unexplained for a long time and the question of the charge transfer in superconductors was left aside till the discovery of the high-$`T_c`$ materials. For these layered materials it was predicted that the superconducting transition induces a charge transfer from CuO2 planes to charge reservoirs. This transfer caused merely by the pairing mechanism has been confirmed by bulk- oriented experiments like the positron annihilation, the $`x`$-ray absorption spectroscopy, and the nuclear magnetic resonance.

Apparently, there are two groups of contradictory experimental results. The pairing contribution is absent in the surface potential but a charge transfer is observed at internal interfaces in the bulk.

As it was indicated recently, there is a charge transfer at the surface which is the interface of superconductor and vacuum. This transfer forms a surface dipole which causes a step $`\phi _\delta `$ in the electrostatic potential. The value of the potential step has been evaluated from the Budd-Vannimenus theorem

$`e\phi _\delta =n{\displaystyle \frac{}{n}}{\displaystyle \frac{f_{\mathrm{el}}}{n}},`$ (4)
where $`f_{\mathrm{el}}`$ is the electronic part of the free energy density (it does not include the electrostatic and magnetic parts). Rickayzen has obtained formula (3) from the general stability condition $`e\phi =\frac{f_{\mathrm{el}}}{n}`$ and the free energy $`f_{\mathrm{el}}=n_s\frac{1}{2}mv^2`$. Using the same free energy one finds that the potential at the surface, $`\phi _0=\phi (0)+\phi _\delta `$, equals the screened Bernoulli potential (2). The surface dipole $`\phi _\delta `$ thus explains the observed magnetic pressure.

In spite of the agreement between observed and theoretically derived voltage one should be reserved about claims that the theory correctly describes the profile of the electrostatic potential in superconductors. Arguments against the theory are similar to those already raised in relation to the interpretation of the measurement of Bok and Klein.

First, the measurement of Morris and Brown explores the entire range of magnetic fields from low up to critical values. The free energy employed by Rickayzen, however, applies in the limit of low magnetic fields only.

Second, materials studied by Morris and Brown are type-I and weak type-II materials so that their behavior is not fully covered by the London theory. As we have shown recently, even for low magnetic fields the electrostatic potential depends on the Ginzburg-Landau (GL) parameter $`\kappa `$. Rickayzen’s formula is recovered for the extreme type-II superconductor $`\kappa \mathrm{}`$. For measured materials with $`\kappa 1/\sqrt{2}`$ the potential $`\phi (0)`$ at the surface is reduced by a factor $`1/3`$ compared to Rickayzen’s formula.

Third, the surface dipole $`\phi _\delta `$ derived from the free energy that covers only the low-field perturbation in the London approximation has the same shortcomings as Rickayzen’s formula. Apparently, $`\phi (0)`$ obtained from Rickayzen’s formula can be far from the correct value and the same applies to $`\phi _\delta `$. Since the sum agrees with the experimental result, one can see that eventual errors tend to compensate each other in the resulting electrostatic potential. In this sense, the surface dipole $`\phi _\delta `$ is consistent with the internal potential $`\phi `$, since both are evaluated using the same free energy.

### D Plan of the paper

As demonstrated for Rickayzen’s theory, the internal electrostatic potential $`\phi `$ and the surface dipole $`\phi _\delta `$ ought to be derived from the same free energy. We have shown that for type-I and weak type-II superconductors, the GL theory yields the internal electrostatic potential which is quite different from the one predicted by Rickayzen’s formula. In this paper we derive the surface dipole within the GL theory.

As in Ref.  we use the Budd-Vannimenus theorem. Here we employ this identity within Bardeen’s extension of the GL theory. Bardeen’s extension offers two advantages. First, it naturally interpolates between the GL theory close to $`T_c`$ and the London theory at low temperatures. Second, it uses material parameters of the Gorter-Casimir two-fluid model, which have a transparent density dependence. In contrast, the parameters of the original GL theory are introduced in the limit $`TT_c`$ and $`T`$ is replaced by $`T_c`$ wherever possible. Since $`T`$ is an independent thermodynamic variable while $`\frac{T_c}{n}0`$, one has to be careful when taking density derivatives. To evaluate the density dependence of the GL parameters, one has to recall the microscopic theory of Gorkov and take the density derivatives of the corresponding parameters before the limit $`TT_c`$ is applied.

The plan of the paper is as follows. In Sec. II we assume temperatures close to $`T_c`$ and an infinitesimally weak magnetic field. In this limit we derive a modification of the Budd-Vannimenus theorem which takes into account a non-zero charge density near the surface. In Sec. III we discuss a general system in the Meissner state and show that the formula for the surface dipole derived in Sec. II applies to any temperature and magnetic field below critical values. Section IV includes a summary.

## II Surface dipole

Let us first estimate the thickness of the surface dipole from thermodynamic considerations. The pairing correlation is weaker on the surface than in the bulk, what results in forces pulling the Cooper pairs inside. Such forces are always balanced by the electrostatic field. The full understanding of this effect will require microscopic studies which are not yet feasible. From the BCS studies it is known, however, that close to the surface on the scale of the BCS coherence length $`\xi _0`$, the gap profile differs from the value given by the GL theory. We thus expect that the surface dipole is somehow linked to this ‘microscopic’ modulation of the gap profile.

To introduce the surface dipole on an intuitive level, let us assume that the system is close to the critical temperature. In this regime, the London penetration depth $`\lambda `$ and the GL coherence length $`\xi `$ are much larger than the BCS coherence length $`\xi _0`$. Since the electrostatic potential induced by the diamagnetic current extends on scales of $`\lambda `$ and $`\xi `$ from the surface, the surface dipole is very narrow on these scales.

We can then define an intermediate scale $`L`$ such that $`\xi _0L\xi ,\lambda `$, as sketched in Fig. 1. On the scale $`L`$, the GL wave function changes only negligibly, i.e., $`\psi (x)\psi (x0)\psi (0)`$ for $`0<x<L`$. We note that the GL boundary condition $`_x\psi =0`$ supports the slow change of $`\psi `$ close to the surface.

As shown by de Gennes, the GL wave function is linearly proportional to the BCS gap, except for the surface region on the scale of the BCS coherence length $`\xi _0`$. Following microscopic theories giving the electrostatic potential in terms of the BCS gap, we expect the electrostatic potential to have similar features, see Fig. 1.

For $`\xi _0x<L`$, the potential is well described by the GL value $`\phi (x)`$. Since $`L\lambda ,\xi `$, the GL prediction of the electrostatic potential changes only negligibly in this region and it is convenient to introduce the extrapolated value $`\phi (x)\phi (x0)\phi (0)`$. The extrapolation of the GL potential towards the surface, $`\phi (0)`$, has to be distinguished from the true surface potential $`\phi _0`$. The difference $`\phi _\delta =\phi _0\phi (0)`$ is caused by the surface dipole we aim to evaluate.

### A Budd-Vannimenus theorem at intermediate scale

Close to the critical temperature, we can take the intermediate region $`xL`$ as a homogeneous ‘bulk’ and follow the idea of Budd and Vannimenus. Let us assume a virtual compression of the crystal lattice such that the background or the lattice charge density is removed from the surface layer of an infinitesimal width $`\delta L`$. The perturbation of the lattice charge density in the infinitesimal layer $`0<x<\delta L`$ is $`\delta \rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}`$. The compression leads to an increase of the charge density in the layer $`\delta L<x<L`$, where $`\delta \rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}\delta L/L`$ is selected to conserve the total charge.

Now we recall the basic idea of the Budd-Vannimenus theorem. The lattice charge enters the jellium model of metals as an external parameter. If one changes this external parameter, the situation corresponds to doing work on the system,

$`\delta 𝒲=S{\displaystyle 𝑑x\delta \rho _{\mathrm{lat}}\frac{f}{\rho _{\mathrm{lat}}}}=S{\displaystyle 𝑑x\delta \rho _{\mathrm{lat}}\phi },`$ (5)
where $`f`$ is the density of the free energy including the electrostatic interaction, and $`S`$ is the sample area. According to the Feynman-Hellmann theorem, the change of the electrostatic potential does not contribute to the work up to the first order in $`\delta \rho _{\mathrm{lat}}`$. Now we can proceed with the algebra. We split the integral into three parts, $`(0,\delta L)`$, $`(\delta L,\xi _0)`$ and $`(\xi _0,L)`$. Since $`\delta L`$ is an infinitesimal displacement, the potential in the layer $`0<x<\delta L`$ can be replaced by the surface value $`\phi _0`$. The surface region $`\delta L<x<\xi _0`$ gives a negligible contribution of the order of $`\xi _0/L`$. In the remaining bulk region $`\xi _0<x<L`$, the electrostatic potential is nearly constant and equals $`\phi (0)`$. The work thus reads

$`\delta 𝒲=S\delta L\rho _{\mathrm{lat}}(\phi (0)\phi _0).`$ (6)
The work increases the free energy $``$ of the system

$`\delta 𝒲=\delta ={\displaystyle \frac{(f_{\mathrm{el}}SL)}{(SL)}}S\delta L=\left(f_{\mathrm{el}}+n{\displaystyle \frac{f_{\mathrm{el}}}{n}}\right)S\delta L,`$ (7)
where $`f_{\mathrm{el}}`$ is the spatial density of the electronic free energy. Note that only the change of the electronic part in the ‘bulk’ is assumed. The surface energy does not change as the surface is merely shifted. The magnetic free energy changes negligibly, because the number of electrons in the layer $`L`$ is not changed by the deformation. Since the condensate velocity changes on the scale $`\lambda L`$, changes of the screening current vanish in the limit $`L/\lambda 0`$.

The first term in (7) results from the reduced volume, $`SLS(L\delta L)`$, and the second one from the corresponding increase of the electron density, $`nn(1+\delta L/L)`$. Equating (6) and (7) we obtain a modification of the Budd-Vannimenus theorem,

$$\rho _{\mathrm{lat}}\phi _\delta \rho _{\mathrm{lat}}\left(\phi _0\phi (0)\right)=f_{\mathrm{el}}n\frac{f_{\mathrm{el}}}{n}.$$
(8)
It describes the step of the potential at the surface due to the surface dipole in terms of the free energy. This relation is the main result of the present paper. In the remaining part we demonstrate how (8) can be used and show that it applies also at low temperatures, where the intermediate region cannot be defined.

Formula (8) differs from the original Budd-Vannimenus (BV) theorem in four points. First, in the original BV theorem one evaluates the total potential step at the surface. Its values are of the order of Volts. Here we evaluate only the change of the potential step which appears as the system goes superconductive. The typical magnitude of this change is of the order of nano-Volts. Second, in the BV theorem the surface potential is related to the potential $`\phi _{\mathrm{}}`$ deep in the bulk. In (8), the extrapolation of the internal potential towards the surface, $`\phi (0)`$, appears instead. Third, in order to cover systems at finite temperatures, we use the free energy instead of the ground-state energy. Fourth, within the original Budd-Vannimenus approach, the electron charge density and the lattice charge density have to be equal because of the charge neutrality. In our approach, the density of electronic charge differs locally from the lattice charge density, $`en\rho _{\mathrm{lat}}`$, due to the charge transfer induced by the magnetic fields.

In the limit of extreme type-II superconductors and weak magnetic fields one can compare formula (8) with the surface dipole evaluated in Ref. . In this limit, the free energy simplifies to $`f_{\mathrm{el}}=f_{\mathrm{}}n_s\frac{1}{2}mv^2`$. The second term provides the current induced changes of the surface dipole derived in Ref. . Due to the bulk free energy $`f_{\mathrm{}}`$ one finds from (8) a finite potential step also in the absence of diamagnetic currents.

For the purpose of the plot in Fig. 1, we have used a small magnitude of the potential step at the surface (the short-dashed line is rather close to the full line). In reality, even the current induced part of the step $`\phi _\delta `$ might achieve a magnitude much larger than the internal potential $`\phi (0)`$. A simple estimate of surface step for the experiment of Morris and Brown ranges from few percent for $`T=0.6T_c`$ to about thirty times of the observed potential at $`T=0.97T_c`$, see Ref. . Since we assume the limit $`TT_c`$, the large values are more appropriate.

### B Application within the GL theory

Now we demonstrate how the relation (8) can be used within the GL theory. To this end we introduce the GL free energy

$$f_{\mathrm{el}}=\frac{1}{2m^{}}\left|\left(i\mathrm{}e^{}𝐀\right)\psi \right|^2+f_{\mathrm{cond}},$$
(9)
where $`\psi `$ is the GL wave function, $`𝐀`$ is the vector potential, $`m^{}=2m`$ and $`e^{}=2e`$ are the mass and the charge of the Cooper pair, and $`f_{\mathrm{cond}}`$ is the free energy of the condensate. It can be either the Gorter-Casimir free energy

$`f_{\mathrm{cond}}={\displaystyle \frac{1}{4}}\gamma T_c^2{\displaystyle \frac{2}{n}}|\psi |^2{\displaystyle \frac{1}{2}}\gamma T^2\sqrt{1{\displaystyle \frac{2}{n}}|\psi |^2}`$ (10)
used by Bardeen or the GL free energy

$`f_{\mathrm{cond}}=\alpha |\psi |^2+{\displaystyle \frac{1}{2}}\beta |\psi |^4.`$ (11)
The GL parameters $`\alpha =\gamma (T_c^2T^2)/2n`$ and $`\beta =\gamma T^2/n^2`$ depend on the temperature $`T`$ and the electron density $`n=n_n+2|\psi |^2`$, where $`n_n`$ is the density of normal electrons. Finally, we add the electromagnetic energy so that the free energy reads

$$f=f_{\mathrm{el}}+\phi (\rho _{\mathrm{lat}}+en)\frac{ϵ_0}{2}E^2+\frac{1}{2\mu _0}B^2,$$
(12)
with the magnetic field $`𝐁=\times 𝐀`$ and the electric field $`𝐄=\phi `$.

Variations of the free energy with respect to its independent variables $`𝐀,\phi ,\psi ,n_n`$ yield the equations of motion in Lagrange’s form

$$\frac{f}{\nu }+\frac{f}{\nu }=0.$$
(13)
For $`\nu =A`$ the variational condition (13) yields the Ampere-Maxwell equation, for $`\nu =\phi `$ the Poisson equation, for $`\nu =\psi `$ the GL equation, and for $`\nu =n_n`$ the condition of zero dissipation,

$$e\phi =\frac{f_{\mathrm{el}}}{n}.$$
(14)
This condition allows one to evaluate the electrostatic potential in the bulk of the superconductor. Of course, one can add any constant to the electrostatic potential.

Formula (14) does not cover the surface dipole on the scale $`\xi _0`$, therefore at the surface it provides the extrapolated bulk value $`\phi (0)`$. We can thus use (14) to rearrange the Budd-Vannimenus theorem (8) as

$$\rho _{\mathrm{lat}}\phi _0=f_{\mathrm{el}}+\phi (0)\left(\rho _{\mathrm{lat}}+en\right).$$
(15)
Now all terms on the right hand side are explicit quantities which one obtains within the GL theory extended by the electrostatic interaction.

### C Convenient approximation

In customary GL treatments, the electrostatic potential and the corresponding charge transfer are omitted. For magnetic properties this approximation works very well, since the relative charge deviation, $`\left(\rho _{\mathrm{lat}}+en\right)/\rho _{\mathrm{lat}}`$, is typically of the order of $`10^{10}`$, leading to comparably small corrections in the GL equation. With the same accuracy one obtains the electronic free energy $`f_{\mathrm{el}}`$. Therefore it is possible to evaluate the surface potential using the approximation

$$\phi _0\frac{f_{\mathrm{el}}}{en_{\mathrm{}}}$$
(16)
which follows from (15) if terms proportional to $`\left(\rho _{\mathrm{lat}}+en\right)/\rho _{\mathrm{lat}}`$ are neglected. By $`n_{\mathrm{}}`$ we have denoted the asymptotic value of $`n`$ deep in the bulk, i.e. the density of pairable electrons, $`\rho _{\mathrm{lat}}=en_{\mathrm{}}`$.

Within approximation (16) one does not have to evaluate the potential profile and the related charge inside the superconductor. This is advantageous, in particular for systems of unknown material parameters $`T_c/n`$ and $`\gamma /n`$.

Turning the argument around, from (16) one can see that the electrostatic potential at the surface cannot be used to measure the material parameters $`T_c/n`$ and $`\gamma /n`$. This fact is already known from the experiment of Morris and Brown.

## III Magnetic pressure

So far we have discussed systems close to the critical temperature, when the validity conditions of the GL theory are well satisfied. In many cases, however, the GL theory is used beyond the limits of its nominal applicability. In these cases the GL coherence length $`\xi `$ and/or the London penetration depth $`\lambda `$ are comparable to, or even shorter than the BCS coherence length $`\xi _0`$ so that the intermediate scale $`L`$ cannot be introduced.

It is possible, however, to follow the original formulation of Budd and Vannimenus and define $`L`$ as the sample thickness, i.e. $`L\xi ,\lambda `$. In this case it is necessary to account for the energy of the magnetic field, since the infinitesimal compression $`\delta L`$ shifts the screening layer inwards into the superconductor.

### A Budd-Vannimenus theorem

For $`L\xi ,\lambda `$, the charge removed from the surface is placed deep in the bulk (the region of the scale of $`\xi ,\lambda `$ gives a negligible contribution) so that the work on the charge reads

$`\delta 𝒲=S\delta L(\phi _{\mathrm{}}\phi _0)\rho _{\mathrm{lat}}.`$ (17)
Compared to the previous treatment we have merely replaced the potential close to the surface by the value deep in the bulk. Similarly, the electronic part of the total free energy changes by

$`\delta _{\mathrm{el}}=\left(f_{\mathrm{el}}^{\mathrm{}}+n{\displaystyle \frac{f_{\mathrm{el}}^{\mathrm{}}}{n}}\right)S\delta L.`$ (18)
Finally, the shift of the screening layer by $`\delta L`$ changes the magnetic energy by an amount $`\delta _B=S\delta LB_0^2/2\mu _0`$, given by the magnetic pressure. Here $`B_0`$ is the value of the magnetic field at the surface.

From $`\delta 𝒲=\delta _{\mathrm{el}}+\delta _B`$ follows

$$\rho _{\mathrm{lat}}(\phi _0\phi _{\mathrm{}})=\frac{B_0^2}{2\mu _0}+f_{\mathrm{el}}^{\mathrm{}}+n_{\mathrm{}}\frac{f_{\mathrm{el}}^{\mathrm{}}}{n}.$$
(19)
As $`L`$ is the thickness of the sample one can use the charge neutrality, $`\rho _{\mathrm{lat}}=en_{\mathrm{}}`$. Since the value of the potential deep in the sample is given by condition (14), using $`e\phi _{\mathrm{}}=\frac{f_{\mathrm{el}}^{\mathrm{}}}{n}`$ from (19) one obtains

$$\phi _0=\frac{B_0^2}{2\mu _0en_{\mathrm{}}}\frac{f_{\mathrm{el}}^{\mathrm{}}}{en_{\mathrm{}}}.$$
(20)
The electrostatic potential observed at the surface is thus given by the magnetic pressure as observed by Morris and Brown.

Note that deriving formula (20) we have not used many assumptions about the system. The condition of zero dissipation (14) is a general thermodynamic relation. The Budd-Vannimenus relation (19), however, is limited to systems with a homogeneous jelly-like background charge. This approximation is acceptable for conventional superconductors, where characteristic scales $`\xi _{\mathrm{BCS}},\xi ,\lambda `$ are much larger than the elementary cell of the crystal. The applicability is questionable for the high-$`T_c`$ materials which due to the layered structure and a short coherence length are far from the jellium model.

As noticed already by Bok and Klein,  there is a simple argument for the formula like (20). If one assumes a slab with magnetic fields $`B_L`$ and $`B_R`$ on the left/right sides, the voltage difference gives

$`\rho _{\mathrm{lat}}(\phi _0^L\phi _0^R)={\displaystyle \frac{1}{2\mu _0}}(B_L^2B_R^2).`$ (21)
The left hand side of this relation represents the electrostatic force (per unit area) on the lattice,

$`F_{\mathrm{elst}}={\displaystyle _L^R}𝑑xE\rho _{\mathrm{lat}}=\rho _{\mathrm{lat}}(\phi _0^L\phi _0^R).`$ (22)
The right hand side is the Lorentz force $`F_{\mathrm{Lor}}=BJ`$ with the mean magnetic field $`B=\frac{1}{2}(B_L+B_R)`$ and the net current $`J=_L^R𝑑xj`$ given by Ampere’s rule, $`B_LB_R=\mu _0J`$. Since the electrostatic field provides the only mechanism by which the force is passed from the electrons to the lattice, the two forces have to be equal, $`F_{\mathrm{Lor}}=F_{\mathrm{elst}}`$. This argument was, however, overlooked in the later studies.

### B Test of the surface relation

The Budd-Vannimenus theorem provides the electrostatic potential (20) in terms of the magnetic field with no regards to the actual potential inside the superconductor. To link formula (20) with the more intuitive derivation from Sec. II, we show that (15) results in the surface potential (20) for any temperature.

For the assumed geometry, the GL equation has an integral of motion, see Bardeen. This integral can be obtained quite generally by the Legendre transformation of the free energy,

$$g=f\underset{\nu }{}\frac{f}{\nu }\nu .$$
(23)
Indeed, if the fields $`\nu `$ obey the equations of motion (13), the gradient $`g=0`$ vanishes, i.e. $`g=\mathrm{const}`$. Deep in the bulk all gradients vanish, therefore $`g=f_{\mathrm{el}}^{\mathrm{}}`$.

From equations (9-12) one finds

$$\underset{\nu }{}\frac{f}{\nu }\nu =\frac{\mathrm{}^2}{m^{}}|\psi |^2ϵ_0E^2+\frac{B^2}{\mu _0}.$$
(24)
With the help of (23-24) and definition (12), one can express the electronic free energy as

$$f_{\mathrm{el}}=f_{\mathrm{el}}^{\mathrm{}}+\frac{\mathrm{}^2}{m^{}}|\psi |^2\frac{ϵ_0E^2}{2}+\frac{B^2}{2\mu _0}\phi (\rho _{\mathrm{lat}}+en).$$
(25)
At the surface, the GL boundary condition demands that $`\psi =0`$ what implies $`E=0`$. The free energy at the surface thus reads

$$f_{\mathrm{el}}=f_{\mathrm{el}}^{\mathrm{}}+\frac{B^2}{2\mu _0}\phi (0)(\rho _{\mathrm{lat}}+en).$$
(26)
From (26) and the surface relation (15) it follows that $`\rho _{\mathrm{lat}}\phi _0=B^2/(2\mu _0)+f_{\mathrm{el}}^{\mathrm{}}`$. This value is identical to (20).

Apparently, we can reverse the procedure. Starting from the general Budd-Vannimenus relation (20) and the general integral of motion (23), we can derive the surface relation (15). Accordingly, the surface relation holds for any temperature, provided that the free energy is a functional of the GL wave function and its first derivative only, $`ff[\psi ,\psi ]`$. This functional can be an arbitrary one.

Perhaps we should explain why we have derived the surface dipole from the Budd-Vannimenus theorem on the intermediate scale, although the more general derivation from the integral of motion is available. There are two reasons. First, the intermediate scale provides at least a qualitative picture of the potential in the vicinity of the surface. This picture might be helpful if measurements sensitive to layers close to the surface will be designed.

Second, within the intermediate scale the surface dipole is treated as a property of the superconducting condensate, what encourages us to hope that formula (15) or its approximation (16) can be used to obtain the surface potential also for cases when the magnetic field has a component perpendicular to the surface. In particular, we expect that it will be applicable to the superconductors in the mixed state, especially to evaluate the electric field generated by vortices penetrating the surface.

## IV Conclusions

In conclusion, the Budd-Vannimenus theorem was modified so that it is applicable to the surface of a superconductor. It allows one to evaluate the electrostatic potential on the surface from the free energy and the bulk electrostatic potential nearby. Formula (16) offers the approximation of the surface potential from the free energy without the actual knowledge of the bulk potential.

For plain surfaces we have recovered the experimentally established fact that the electrostatic potential equals the magnetic pressure divided by the density of pairable electrons. This experimental law was confirmed also for type-I and weak type-II superconductors, while the previous theoretical treatments were restricted to weak magnetic fields and extreme type-II superconductors. The presented theory is free of these limitations.

It was shown that thermodynamic corrections do not influence the surface electrostatic potential, measurable e.g. via contactless capacitive pickup. Consequently, contrary to earlier expectations, the density dependence of the critical temperature cannot be estimated in this way. On the other hand, the relation between the surface electrostatic potential and the magnetic pressure shows, that such a measurement allows one to determine the density of charge carriers without knowledge of any other material parameters.

In this paper we have derived only the amplitude of the potential step. The detailed profile of the electrostatic potential including its modulation at the surface can be obtained by a microscopic approach like the Bogoljubov-de Gennes theory extended recently to cover the electrostatic phenomena. For microscopic calculations, the Budd-Vannimenus theorem can serve as a test of accuracy of the numerical procedure, similarly as it is used in the theory of metal surfaces.

###### Acknowledgements.

This work was supported by MŠMT program Kontakt ME601 and GAČR 202/03/0410, GAAV A1010312 grants. The European (ESF) program VORTEX is also acknowledged.