document: there are many problems concerning electronic structure where attention is focussed on a small region of a larger system , at surfaces or defects in crystals being perhaps the most common . 
   let us call this region i , figure [ fig:1 ] , and the rest of the system region ii . 
   although not of primary interest region ii can not be ignored , since in general the electron wave functions in i will be sensitive to the contents of region ii . 
   some time ago inglesfield @xcite derived an embedding scheme which enables the single - particle schrdinger equation to be solved explicitly only in region i. the influence of region ii is taken into account exactly by adding an energy - dependent non - local potential to the hamiltonian for region i , which constrains the solutions in i to match onto solutions in ii . 
   this embedding method has been developed into a powerful tool most notably for surface electronic structure problems @xcite where it has found widespread application especially to situations where an accurate description of the spectrum of electron states is necessary . 
   examples include studies of image states @xcite , surface states at metals surfaces @xcite , static and dynamic screening @xcite , atomic adsorption and scattering at surfaces @xcite , studies of surface optical response @xcite and field emission @xcite . 
   recent applications to transport problems have also been described @xcite . for a review of the embedding method 
   see inglesfield @xcite . 
   = 80 mm    in the case of materials containing heavier elements , relativistic effects can be significant @xcite and lead to important deviations from the electronic structure as predicted by the schrdinger equation  shifts in inner core levels of 5d elements are typically several 100 or 1000 ev , valence bands shifts are on the ev scale and spin orbit splitting is often measured in tenths of ev . even ignoring the concomitant changes in electron wave functions these shifts can reorder levels and so affect calculated densities , fundamental to the determination of ground state properties within the density functional framework @xcite . 
   for this reason most of the conventional electronic structure techniques developed for accurately solving the single - particle schrdinger equation in solids have subsequently been modified to deal with the dirac equation , including the relativistic augmented plane wave method @xcite , relativistic linear muffin - tin orbital method @xcite , relativistic augmented spherical wave method @xcite and the relativistic multiple - scattering method @xcite , and each has subsequently been used in studying a diverse range of problems . 
   the last method alone has formed the basis of calculations of photoemission @xcite , magnetocrystalline anisotropy @xcite , hyperfine interactions @xcite and magnetotransport @xcite amongst other topics . 
   inglesfield s embedding method has particular advantages that encourage its extension to the relativistic case . 
   it permits the inclusion of extended substrates for surface and interface calculations , enables the study of isolated point defects in solids and being a basis set technique is highly flexible and permits full - potential studies with relative ease . at surfaces extended substrates ( as against the use of the supercell or thin - film approximation in which the crystal is approximated by a small number of layers , typically 5 - 7 ) enable the proper distinction between surface states , resonances and the continuum of bulk states @xcite . 
   the behaviour of the w(110 ) surface @xcite where the addition of half a monolayer of li is observed to _ increase _ the spin - orbit splitting of a surface state by @xmath2 ev ( resulting in fermi surface crossings separated by @xmath3% of the brillouin zone dimension ) typifies a type of problem a relativistic embedding scheme could address . 
   indeed each of the topics mentioned at the end of the previous paragraph are relevant at surfaces and/or interfaces , and could be usefully investigated within a relativistic embedding framework .    in this paper 
   we develop an embedding scheme for the dirac equation that parallels inglesfield s scheme for the schrdinger equation . 
   inglesfield s starting point is the expectation value of the hamiltonian using a trial wave function which is continuous in amplitude but discontinuous in derivative across the surface @xmath1 separating i and ii . 
   the first order nature of the dirac equation precludes the use of a similar trial function . 
   instead , in the following section we use a trial function in which the large component is continuous and the small component discontinuous across @xmath1 . 
   continuity in the small component is restored when the resulting equations are solved exactly . using the green function for region 
   ii we are able to derive an expression for the expectation value purely in terms of the trial function in i. in section [ section : app ] the application of the method is illustrated by calculating the eigenstates of a hydrogen atom within a cavity and in section [ section : green ] we determine the green function for the embedded region . 
   section [ section : monolayer ] briefly illustrates the method applied to a sandwich structure where relativistic effects are marked . 
   we conclude with a brief summary and discussion . 
   in this section we consider region i joined onto region ii ( figure [ fig:1 ] ) , and derive a variational principle for a trial wave function @xmath4 defined explicitly only within region i. we are primarily interested in the positive energy solutions of the dirac equation @xcite , and so we refer to the upper and lower spinors of the dirac bi - spinor solutions as the large and small components of the wave function respectively . 
   we notionally extend @xmath4 into ii as @xmath5 , an exact solution of the dirac equation at some energy @xmath6 , with the large components of @xmath4 and @xmath5 ( @xmath7 and @xmath8 ) matching on the surface @xmath1 separating i and ii , but with no constraint upon the small components ( @xmath9 and @xmath10 ) , figure [ fig:2 ] . 
   the expectation value for the energy @xmath11 is then @xmath12 } { \langle\varphi|\varphi\rangle_{\mathrm{i}}+ \langle\chi|\chi\rangle_{\mathrm{ii } } } \label{eqn : exp}\ ] ] where @xmath13 . 
   ( for clarity we omit the interaction @xmath14 which appears in the relativistic density functional theory @xcite neglecting orbital and displacement currents , where @xmath15 is a `` spin - only '' effective magnetic field containing an external and exchange - correlation contribution . 
   its inclusion has no consequences for the derivation . ) 
   the first two terms in the numerator are the expectation value of the hamiltonian through regions i and ii , and the third the contribution due to the discontinuity in the small component of the wave function on @xmath1 ( in this and the following , surface normals are directed from i to ii ) . 
   = 60 mm    we eliminate reference to @xmath5 by introducing two relations . 
   firstly , for @xmath16 , @xmath5 satisfies the dirac equation at energy @xmath6 @xmath17\chi=0 \label{eqn : dirii}\ ] ] and differentiating with respect to @xmath6 the energy derivative of @xmath5 , @xmath18 , satisfies @xmath19\dot{\chi}=\chi\qquad \bi{r}\in\mathrm{ii}. \label{eqn : eder}\ ] ] multiplying the hermitian conjugate of the first equation by @xmath20 from the right , multiplying the second from the left by @xmath21 , subtracting and integrating over region ii gives a relation between the normalisation of @xmath5 in ii and the amplitude on @xmath1 : @xmath22 we have assumed that @xmath5 vanishes sufficiently strongly at infinity .    for the second relation 
   we introduce the green function ( resolvant ) @xmath23 corresponding to equation ( [ eqn : dirii ] ) : @xmath24g=-\delta(\bi{r}-\bi{r } ' )   \qquad \bi{r},\bi{r}'\in \mathrm{ii}. \label{eqn : gf}\ ] ] multiplying the hermitian conjugate of this equation by @xmath5 from the right , and subtracting @xmath25 times equation ( [ eqn : dirii ] ) , integrating over region ii and then using the reciprocity of the green function gives @xmath26 we see that the green function relates the amplitude of the wave function on @xmath1 to the amplitude at any point within ii . 
   in particular , we can obtain a relation between the large and small components of @xmath5 on s. writing the @xmath27 green function as @xmath28 where each entry is a @xmath29 matrix , substituting into equation ( [ eqn : chis ] ) , and rearranging the two equations coupling the small and large components of @xmath5 gives @xmath30 where @xmath31 it follows from ( [ eqn : chis ] ) that the green functions in ( [ eqn : gam ] ) are the limiting forms of @xmath32 as @xmath33 from within ii . 
   equations ( [ eqn : norm ] ) and ( [ eqn : cgc ] ) are the desired results that enable us to express the expectation value @xmath11 in ( [ eqn : exp ] ) in terms of @xmath4 alone . 
   after substitution and use of the continuity of the large components @xmath34 on @xmath1 we obtain @xmath35 } { \langle\varphi|\varphi\rangle_{\mathrm{i}}- c^2\hbar^2 \int_s \rmd\bi{r}_s\cdot \varphi_{\rm l}^\dag \bsigma \int_s \rmd\bi{r}'_s \cdot \dot{\gamma } \bsigma\ , \varphi_{\rm l } } . 
   \label{eqn : exp2}\ ] ] this is an expression for the expectation value of the energy @xmath11 , given purely in terms of the trial function @xmath4 in region i and on the surface @xmath1 , with all details of region ii entering via @xmath36 and its energy derivative . following the convention in the non - relativistic embedding scheme we shall refer to @xmath36 as the embedding potential .    to see what this variational principle means in practice , we consider variations in @xmath37 , whereby @xmath38 } { \langle\varphi|\varphi\rangle_{\mathrm{i}}- c^2\hbar^2\int_s \rmd\bi{r}_s\cdot \varphi_{\rm l}^\dag \bsigma \int_s \rmd\bi{r}'_s \cdot \dot{\gamma } \bsigma\ , \varphi_{\rm l } } , \nonumber\\ \label{eqn : de1}\end{aligned}\ ] ] so that solutions @xmath4 stationary with respect to arbitrary variations @xmath39 satisfy @xmath40 the first expression indicates @xmath4 is a solution of the dirac equation at energy @xmath11 in region i. comparing the second with ( [ eqn : cgc ] ) shows that @xmath4 also possesses the correct relationship between large and small components on @xmath1 , the surface separating i and ii , to match onto solutions in ii . 
   the term @xmath41 provides a first order correction to @xmath42 so that the boundary condition is appropriate for energy @xmath11 .    in practice expression ( [ eqn : exp2 ] ) may be used to obtain solutions of the dirac hamiltonian by inserting a suitably parameterised trial function and varying the parameters to obtain a stationary solution . 
   this is conveniently achieved by expanding the trial solution in a finite basis of separate large and small component spinors @xmath43 + \sum_{n=1}^{n_{\rm s } } a_{{\rm s},n }   \left[\begin{array}{c}0\\ \psi_{{\rm s},n}(\bi{r})\end{array}\right ] = \left[\begin{array}{cc } \bpsi_{\rm l}(\bi{r } ) & 0 \\ 0 & \bpsi_{\rm s}(\bi{r } ) \end{array}\right ] \left [ \begin{array}{c } \bi{a}_{\rm l } \\ \bi{a}_{\rm s } \end{array } \right].\ ] ] the matrix in the final expression is @xmath44 by @xmath45 , and the column vector contains the @xmath45 coefficients . 
   substituting into ( [ eqn : exp2 ] ) we find states @xmath46 that are stationary with respect to variations in the expansion coefficients @xmath47 are then given by the eigenstates of a generalised eigenvalue problem of the form @xmath48 \left [ \begin{array}{c } \bi{a}_{\rm l } \\ 
   \bi{a}_{\rm s } \end{array } \right ] = w \left [ \begin{array}{cc } { o}_{\rm ll } & 0 \\ 0 & { o}_{\rm ss } \end{array } \right ] \left [ \begin{array}{c } \bi{a}_{\rm l } \\ 
   \bi{a}_{\rm s } \end{array } \right ] , \label{eqn : a1}\ ] ] where @xmath49_{nn'}= \int_i \psi_{{\rm l},n}^\dag(\bi{r } ) \left(v(\bi{r})+mc^2 \right )   \psi_{{\rm l},n'}(\bi{r } ) \rmd\bi{r } \nonumber \\ 
   \fl \phantom{\left[{h}_{\rm ll}\right]_{nn'}=}{}+c^2\hbar^2\!\int_s \rmd\bi{r}_s\cdot   \psi_{{\rm l},n}^\dag(\bi{r}_s)\bsigma\ ! \int_s \rmd\bi{r}'_s \!\cdot\ !   \left [ \gamma(\bi{r}_s,\bi{r}'_s;w)\ ! - \ ! 
   w \dot{\gamma}(\bi{r}_s,\bi{r}'_s;w)\right ] \bsigma\psi_{{\rm l},n'}(\bi{r}'_s ) \label{eqn : h2}\\ \fl \left[{h}_{\rm ls}\right]_{nn'}= \int_i \psi_{{\rm l},n}^\dag(\bi{r } ) c\bsigma\cdot \widehat{\bi{p } }   \psi_{{\rm s},n'}(\bi{r } ) \rmd\bi{r}+\rmi c\hbar\int_s \rmd\bi{r}_s\cdot   \psi_{{\rm l},n}^\dag(\bi{r}_s)\bsigma\psi_{{\rm s},n'}(\bi{r}_s)\label{eqn : h3}\\ \fl \left[{h}_{\rm sl}\right]_{nn'}= \int_i \psi_{{\rm s},n}^\dag(\bi{r } ) c\bsigma\cdot\widehat{\bi{p } } \psi_{{\rm l},n'}(\bi{r } ) \rmd\bi{r } \label{eqn : h4}\\ \fl \left[{h}_{\rm ss}\right]_{nn'}= \int_i \psi_{{\rm s},n}^\dag(\bi{r } ) \left(v(\bi{r})-mc^2 \right )   \psi_{{\rm s},n'}(\bi{r } ) \rmd\bi{r } \label{eqn : h5}\\ \fl \left[{o}_{\rm ll}\right]_{nn'}= \int_i \psi_{{\rm l},n}^\dag(\bi{r } )   \psi_{{\rm s},n'}(\bi{r } ) \rmd\bi{r}\nonumber\\ \fl \phantom{\left[{o}_{\rm ll}\right]_{nn'}= } { } -c^2\hbar^2\int_s \rmd\bi{r}_s\cdot   \psi_{{\rm l},n}^\dag(\bi{r}_s)\bsigma \int_s \rmd\bi{r}'_s \cdot   \dot{\gamma}(\bi{r}_s,\bi{r}'_s;w ) \bsigma\psi_{{\rm l},n'}(\bi{r}'_s ) \label{eqn : h6}\\ \fl \left[{o}_{\rm ss}\right]_{nn'}= \int_i \psi_{{\rm s},n}^\dag(\bi{r } ) \psi_{{\rm s},n'}(\bi{r } )   \rmd\bi{r}.\label{eqn : h7}\end{aligned}\ ] ] of course the spectrum of the dirac hamiltonian is unbounded below , and care must be taken to prevent solutions collapsing to negative energies . 
   this can be avoided through the use of a kinetically balanced basis @xcite in which there is a one - to - one relationship between large and small component spinors , @xmath50 , and where the small component spinors are given by @xmath51 the upper half of the spectrum of the @xmath52 eigenstates of ( [ eqn : a1 ] ) then provide approximations to the spectrum of electronic states . 
   = 80 mm    to illustrate the application of the relativistic embedding scheme we consider a model problem of a hydrogen atom within a spherical cavity , finding bound states of the dirac equation corresponding to the potential illustrated in figure [ fig:3 ] : @xmath53 where @xmath54 and @xmath55 . 
   we choose this model as the bound states may also be found straightforwardly by alternative methods . 
   region i , the region to be treated explicitly , is the sphere of radius @xmath56 centered on @xmath57 . 
   the external region ii where @xmath58 is replaced by an embedding potential acting on the surface of the sphere . 
   the value of the embedding potential is most readily evaluated from equation ( [ eqn : cgc ] ) . a general solution to the dirac equation at some energy @xmath6 in region 
   ii and satisfying the appropriate boundary conditions is @xcite @xmath59 where @xmath60 , @xmath61 , @xmath62 a spin - angular function , @xmath63 a modified spherical bessel function of the third kind @xcite , @xmath64 , @xmath65 and @xmath66 the spherical symmetry of region ii means the the embedding potential @xmath36 may be expanded on @xmath1 as @xmath67 and substituting ( [ eqn : a2b ] ) and ( [ eqn : a2 ] ) into ( [ eqn : cgc ] ) leads to @xmath68 using ( [ eqn : gam ] ) with the green function for constant potential ( @xmath69 ) @xmath70 gives the same result but after rather more involved manipulations . 
   @xmath71 is a modified spherical bessel function of the first kind .    because of the spherical symmetry we can determine separately states with a given angular character @xmath72 . using as a basis set for the large component spinors @xmath73 so that the small component spinors ensuring kinetic balance are @xmath74 e^{-r } , \label{eqn : bas2}\ ] ] the matrix elements become @xmath75_{nn'}= \int_0^r g_{n}(r)\left[-\frac{\lambda}{r}+mc^2\right ] g_{n'}(r ) \rmd r 
   \nonumber\\ + \hbar^2c^2r^2g_{n}(r)g_{n'}(r ) \left[\gamma_\kappa(w)-w\dot{\gamma}_\kappa(w)\right ] \\ 
   \fl \left[{h}^{(\lambda)}_{\rm ss}\right]_{nn'}= \int_0^rf_{n\kappa}(r)\left[-\frac{\lambda}{r}-mc^2\right]f_{n'\kappa}(r )   \rmd r\\ \fl \left[{h}^{(\lambda)}_{\rm ls}\right]_{nn'}=-\hbar c 
   \int_0^r g_{n}(r)\left[\frac{\rmd f_{n'\kappa}(r)}{\rmd r}- \frac{\kappa}{r}f_{n'\kappa}(r)\right ] \rmd r + \hbar c g_{n}(r)f_{n'\kappa}(r)\\ \fl \left[{h}^{(\lambda)}_{\rm sl}\right]_{nn'}=\hbar c 
   \int_0^r f_{n\kappa}(r)\left[\frac{\rmd g_{n'}(r)}{\rmd r}+ \frac{\kappa}{r}g_{n'}(r)\right ] \rmd r \\ 
   \fl \left[{o}^{(\lambda)}_{\rm ll}\right]_{nn'}= \int_0^r g_n(r)g_{n'}(r ) \rmd r -\hbar^2c^2r^2g_n(r)g_{n'}(r)\dot{\gamma}_\kappa(w ) 
   \\ \fl \left[{o}^{(\lambda)}_{\rm ss}\right]_{nn'}= \int_0^r f_{n\kappa}(r)f_{n'\kappa}(r ) \rmd r\end{aligned}\ ] ] the eigenvalues only depend upon the quantum number @xmath76 . in table 
   [ table:1 ] the lowest two eigenvalues of @xmath77 symmetry ( corresponding to the @xmath78 and @xmath79 of free hydrogen ) are shown as a function of basis set size and for different values of the energy @xmath6 at which the embedding potential is evaluated , for the case @xmath80 , @xmath81 . 
   for comparison also given are the values found by matching the external solution ( [ eqn : a2 ] ) to the regular internal solution , which can be expressed in terms of confluent hypergeometric functions@xcite . for 
   a given fixed @xmath6 the eigenvalues converge from above to values that are equal or above the exact values . 
   the further @xmath6 lies from the eigenvalue , the larger the difference between the limiting value for large basis sets and the correct value . 
   however , the influence of the @xmath82 terms in ( [ eqn : exp2 ] ) means the error is relatively small . when @xmath83 , the lowest eigenvalue found with @xmath84 is @xmath85 ha and in error by only 0.0000044 ha , a factor @xmath86 smaller than the error in @xmath6 . 
   @cccc @xmath87 & @xmath88 & @xmath83 & @xmath89 + 2 & -0.4111620 ,  1.6980995 & -0.4111527 , 1.6949300 & -0.4111624 , 1.6896482 + 4 & -0.4451482 ,  0.9129418 & -0.4451439 , 0.9126817 & -0.4396204 , 0.9714775 + 6 & -0.4455519 ,  0.8914789 & -0.4455477 , 0.8912219 & -0.4455520 , 0.8910268 + 8 & -0.4455532 ,  0.8912708 & -0.4455488 , 0.8910141 & -0.4455532 , 0.8908194 
   + `` exact '' & -0.4455532 , 0.8908194 & -0.4455532 , 0.8908194 & -0.4455532 , 0.8908194 +    differentiating ( [ eqn : exp2 ] ) with respect to the trial energy @xmath6 shows the expectation value is stationary at @xmath90 . in this case @xmath11 is given by the solutions of @xmath91 } { \langle\varphi|\varphi\rangle_{\mathrm{i } } } . 
   \label{eqn : exp3}\ ] ] eigenfunctions @xmath4 solving this equation satisfy the dirac equation within i and the relationship between small and large components on @xmath1 ( [ eqn : de2 ] ) is exact . the final column in table [ table:1 ] 
   shows the lowest two positive energy eigenvalues of ( [ eqn : exp3 ] ) , again as a function of basis set size . 
   the eigenvalues again converge from above , and by @xmath84 reproduce the exact values by at least 7 significant figures . 
   it is worth noting that with this particular basis set increasing @xmath87 much further leads to some numerical difficulties due to overcompleteness . for more accurate 
   work a more suitable basis set should be used . 
   it should also be noted that conventional finite basis set calculations using a basis satisfying kinetic balance can given eigenvalues that lie below exact limiting values by an amount of order @xmath92 @xcite , and similar behaviour is expected in this embedding scheme . 
   most practical applications of the schrdinger embedding scheme have actually used the green function of the embedded system . 
   this is a more convenient quantity when dealing with systems where the spectrum is continuous , such as at surfaces or defects in solids . 
   we therefore consider the green function for the embedded dirac system 
   .    differentiating ( [ eqn : exp2 ] ) with respect to @xmath6 shows @xmath11 is stationary when @xmath90 , as would be expected . in this case 
   stationary solutions satisfy the embedded dirac equation @xmath93 where , introducing @xmath94 , the component of @xmath95 in the direction normal to the surface @xmath1 ( from i to ii ) at @xmath96 , the additional term @xmath97 enforcing the embedding is @xmath98 . 
   \label{eqn : gf2}\end{aligned}\ ] ] the corresponding green function satisfies @xmath99 for @xmath100 . 
   a similar line of argument to that given by inglesfield @xcite for the embedded schrdinger equation shows that this green function is identical for @xmath100 to the green functions @xmath101 for the entire system i@xmath102ii . for simplicity assuming i@xmath102ii constitute a finite system so that the spectrum is discrete , the green function @xmath101 is given by @xmath103 where @xmath104 is the eigenvalue corresponding to eigenstate @xmath105 of the entire system , normalised to unity over i@xmath102ii . for a given @xmath11 , the green function solving ( [ eqn : gf3 ] ) 
   can be expanded in terms of the eigenstates @xmath106 of the corresponding homogeneous equation @xmath107 normalised to unity over i , as @xmath108 clearly @xmath109 has poles at @xmath110 . at these energies ( [ eqn : gf5 ] ) becomes the exact embedded dirac equation ( [ eqn : gf1 ] ) so as we have seen the poles will occur at eigenstates of the entire system and the spectrum of @xmath109 and @xmath101 coincide . it remains to show the poles of @xmath109 have the appropriate weight . 
   the residue of @xmath109 at @xmath111 is @xmath112 the second term in the denominator is precisely the additional factor necessary to correctly normalise the states ( see ( [ eqn : norm ] ) , ( [ eqn : cgc ] ) ) so that @xmath113 the residues of the green function of the embedded system and those of the entire system are identical . 
   hence the two green functions are identical for @xmath114 i.    for practical calculations the green function can be expanded using a double - basis of separate large and small component spinors : @xmath115{g}(w ) \left[\begin{array}{cc } \bpsi_{\rm l}(\bi{r } ' ) & 0 \\ 0 & \bpsi_{\rm s}(\bi{r } ' ) \end{array}\right]^\dag.\ ] ] the matrix elements of the matrix of coefficients @xmath116 may be found by substituting into ( [ eqn : gf3 ] ) , multiplying from the right by the vector of basis functions , multiplying from the left by the hermitian transpose of the vector of basis functions , and integrating over region @xmath117 . 
   this leads to @xmath118^{-1}\ ] ] where the overlap and hamiltonian matrices have their previous definitions ( [ eqn : h2]-[eqn : h7 ] ) with @xmath119 .    as an illustration 
   we calculate the local density of states for the confined hydrogen model at energies above @xmath120 where the spectrum is continuous . integrating over the embedded region this 
   is given by @xmath121 figure [ fig:4 ] shows the @xmath122 wave local density of states for @xmath80 , @xmath123 , calculated with varying number of basis functions . 
   the basis functions ( [ eqn : bas1 ] ) , ( [ eqn : bas2 ] ) are not particularly appropriate for representing the continuum wave solutions , and so convergence is only achieved using a relatively large set ; however the results serve to illustrate the systematic improvement that accompanies an increasing number of basis functions . 
   the local density of states shows two resonances , the precursors of bound states that exist when any of @xmath56 , @xmath124 or @xmath120 are increased sufficiently . 
   as a further example , one that provides a test of the relativistic embedding scheme when applied to a more challenging problem , we use it to calculate the local density of states on a silver monolayer in a au(001)/ag / au(001 ) sandwich structure . using the embedding scheme only the region occupied by the ag monolayer 
   is explicitly treated . 
   this is region i , with the two au halfspaces to either side entering the calculation via embedding potentials expanded on planar surfaces . 
   then , using bloch s theorem the calculation is performed within a unit cell containing one atom . 
   the full technical details will be described elsewhere , but briefly the green function at two - dimensional wave vector @xmath125 is expanded in a set of linearised augmented relativistic plane waves . 
   we use large component basis functions @xmath126\omega_\lambda(\widehat{\bi{r } } ) & & \bi{r } \in \mbox{muffin - tin } \end{array } \right.\ ] ] where @xmath127 is a pauli spinor , @xmath128 , with @xmath129 a two - dimensional reciprocal lattice vector and @xmath130 , @xmath131 and where @xmath132 exceeds the width of the embedded region ensuring variational freedom in the basis . 
   the function @xmath133 is the large component of the wavefunction that satisfies the radial dirac equation for the spherically symmetric component of the potential at some pivot energy ; @xmath134 is the energy derivative of @xmath135 . 
   the matching coefficients @xmath136 , @xmath137 ensure continuity of the basis function in amplitude and derivative at the muffin - tin radius . 
   the small component basis functions are chosen to satisfy kinetic balance . 
   overlap and hamiltonian matrix elements follow directly from these basis functions 
   . the embedding potential is obtained from ( [ eqn : cgc ] ) using the general expression for a wavefunction outside a surface at wave vector @xmath125 . 
   this gives for the embedding potential describing the left au half space @xmath138_{\bi{g}\sigma\bi{g}'\sigma'}\nonumber\\ \times \exp(\rmi(\bi{k}+\bi{g})\cdot \bi{r}_s ) \exp(-\rmi(\bi{k}+\bi{g}')\cdot \bi{r}'_s ) \,\varphi_\sigma 
   \otimes \varphi_{\sigma ' } \label{eqn : embpot}\end{aligned}\ ] ] with @xmath139 the reflection matrix @xmath140 is found using standard layer - scattering methods @xcite 
   . a similar approach may be used to obtain an embedding potential for the right half space , which unlike the non - relativistic case differs from that for the left half space . 
   figure [ fig:5 ] compares the local density of states calculated using the relativistic embedding technique for an embedded ag monolayer using embedding potentials corresponding to au(001 ) with that found for an au(001)/ag / au(001 ) sandwich geometry using relativistic scattering theory @xcite . 
   the same au and ag potentials has been used in each case , and the local density of states found within the same muffin - tin volume . 
   therefore the results obtained with the two methods should be comparable , and we find that they are indistinguishable . 
   this confirms that the embedding potential ( [ eqn : embpot ] ) imposes the correct variational constraint upon wave functions for the embedded ag monolayer so that they replicate the behaviour of an extended au(001)/ag / au(001 ) sandwich structure . 
   the inset in figure [ fig:5 ] shows the local density of states in the non - relativistic limit ( @xmath141 ) , indicating the significant relativistic effects on the electronic structure which are correctly reproduced with this dirac - embedding scheme . 
   we have outline above an embedding scheme for the dirac equation . 
   it enables the dirac equation to be solved within a limited region i when this region forms part of a larger system , i@xmath102ii . 
   region ii is replaced by an additional term added to the hamiltonian for region i , and which acts on the surface @xmath1 separating i and ii . 
   the embedding scheme is derived using a trial function in which continuity in the small component across @xmath1 is imposed variationally . 
   expanding the wave function in a basis set of separate large and small component spinors , the problem of variational collapse is avoided by using a basis satisfying kinetic balance . calculating the spectrum of a confined hydrogen atom , the method is shown to be stable and converge to the exact eigenvalues . 
   we have also derived the green function for the embedded hamiltonian and illustrated its use in the continuum regime of the same confined hydrogen system and an au / ag / au sandwich structure . 
   these are demonstration calculations 
   future applications are likely to be to defects and surfaces of materials containing heavier ( typically 5@xmath142 ) elements , within the framework of density functional theory . 
   it is worthwhile to discuss further the use of a trial function that is discontinuous in the small component , since such a wave function gives rise to a discontinuous probability density and so would normally be dismissed in quantum theory . in non - relativistic quantum mechanics 
   discontinuous trial functions are not permitted , since they possess infinite energy . 
   however the dirac equation is first order in @xmath143 and as we have seen a perfectly regular expectation value of @xmath0 results . exploiting this freedom 
   , the embedding scheme outlined above leads to solutions that are continuous in both large and small component _ only _ when the embedding potential @xmath144 is evaluated at the same energy @xmath6 as the energy @xmath11 that appears in the dirac equation itself , for then the relationship between small and large components on @xmath1 inside ( equation [ eqn : de2 ] ) and outside ( equation [ eqn : cgc ] ) coincide , the large components matching by construction . 
   this may be achieved for example via the iterative scheme used in connection with equation ( [ eqn : exp3 ] ) and the final column of table [ table:1 ] , or explicitly when determining the green function as in section [ section : green ] . 
   these are the methods in which the non - relativistic embedding scheme has been most widely used . 
   when @xmath6 and @xmath11 do not coincide , the solutions obtained via this embedding scheme will retain small components that are discontinuous across @xmath1 . 
   this may be unacceptable for certain applications , but the solutions continue to be valid approximations at least in as much as they provide estimates of the energies of the solutions of the dirac equation , and so could suffice e.g. for interpreting spectroscopic measurements . 
   this embedding scheme places no greater emphasis on a discontinuity in the small amplitude at @xmath1 than on an incorrect ( but continuous ) amplitude elsewhere within the embedded region . 
   it aims merely to optimise the energy of the state , and will retain a discontinuity in the small component if in doing so it can better ( in terms of energy ) approximate the solution inside the embedded region . in the non - relativistic embedding scheme 
   the discontinuous derivative of the trial function implies a probability current ( and electric current ) that is discontinuous across the embedding surface . 
   this is similarly unphysical , yet numerous applications such as those cited above have demonstrated the utility and accuracy of the method . indeed , there have been many applications in which this scheme has been used to determine currents and or transport properties , such as in relation to surface optical response @xcite and electron transport in electron waveguides or through domain walls @xcite . 
   the reason for the success of these calculations is that they employed schemes in which the embedding potential was evaluated at the correct energy , ensuring that the derivative of the wavefunction was continuous across @xmath1 . in practise 
   there have been few calculations using the non - relativistic embedding scheme in which the energies did not coincide . 
   there are a number of aspects of the method which are worthy of further consideration . 
   we started with a trial function in which by construction the large component was continuous and the small component discontinuous across the surface @xmath1 dividing i and ii . 
   we could have reversed these conditions , leading to a similar embedded dirac equation but with a modified embedding term . 
   the particular choice was motivated by the wish to have a theory which behaves reasonably in the limit @xmath141 when the small component becomes negligible  a discontinuous amplitude is not permissible in trial solutions to the schrdinger equation . 
   however , the behaviour of the alternative formulation should be investigated . perhaps in connection with this there is the question of the spectrum of negative energy solutions , to which we have paid scant attention . 
   exploring the @xmath141 limit it might be possible to identify how to embed a relativistic region i within a region ii treated non - relativistically 
   a 5@xmath142 overlayer on a simple metal substrate might be a physical system where such a treatment is appropriate 
   . there could be benefits in terms of computational resources expended if the embedding potential could be determined within the framework of a non - relativistic calculation , and there might also be useful insights in terms of simple models . 
   finally , in terms of implementation for realistic systems , some of the novel schemes for deriving embedding potentials @xcite could certainly be adapted to the relativistic case . 
   it would also be worthwhile to consider whether it is possible to use a restricted electron - like basis , in which the large and small component spinors are combined . 
   this is common practice in most relativistic electronic structure calculations for solids when using basis set techniques ( e.g. @xcite ) , and would result in significant computational efficiencies .
summary: an embedding scheme is developed for the dirac hamiltonian @xmath0 . dividing space into regions i and ii separated by surface @xmath1 , 
   an expression is derived for the expectation value of @xmath0 which makes explicit reference to a trial function defined in i alone , with all details of region ii replaced by an effective potential acting on @xmath1 and which is related to the green function of region ii . 
   stationary solutions provide approximations to the eigenstates of @xmath0 within i. the green function for the embedded hamiltonian is equal to the green function for the entire system in region i. application of the method is illustrated for the problem of a hydrogen atom in a spherical cavity and an au(001)/ag / au(001 ) sandwich structure using basis sets that satisfy kinetic balance .