text:
# Introduction

## Introduction

The coupling of the vector and pseudoscalar heavy-light mesons to the pion, $`g_{H^{}H\pi }`$ , and the related $`\widehat{g}_Q`$ defined through111The subscript Q $`\{c,b\}`$ refers to the heavy quark constituent of the heavy-light meson.

$`g_{H^{}H\pi }={\displaystyle \frac{2\sqrt{m_Hm_H^{}}}{f_\pi }}\widehat{g}_Q,`$ (1)
have been extensively studied in the literature, both for phenomenological and theoretical reasons (a list of results can be found in refs. ). $`\widehat{g}_Q`$ is particularly important in the studies of the $`q^2`$-dependence of the form factor $`F_+(q^2)`$ which encodes the non-perturbative QCD dynamics of $`D\pi `$ and $`B\pi `$ semileptonic decays. To better illustrate this, we consider the dispersion formula for $`F_+^{B\pi }(q^2)`$ :

$`F_+^{B\pi }(q^2)={\displaystyle \frac{1}{2m_B^{}}}{\displaystyle \frac{f_B^{}g_{B^{}B\pi }}{1{\displaystyle \frac{q^2}{m_B^{}^2}}}}+{\displaystyle \frac{1}{\pi }}{\displaystyle _{(m_\pi +m_B)^2}^{\mathrm{}}}𝑑t{\displaystyle \frac{\mathrm{Im}F_+(t)}{tq^2i\epsilon }},`$ (2)
where $`m_B^{}`$ and $`f_B^{}`$ are the mass and the decay constant of the $`B^{}`$-meson. Near the $`B^{}`$-pole, $`F_+`$ can be approximated by the first term of the r.h.s of eq. (2), which is why $`g_{B^{}B\pi }`$ is so important. Inversely, from the $`q^2`$-dependence deduced from the lattice data for $`F_+^{B\pi }`$, one can extract indirectly the value of $`g_{B^{}B\pi }`$ (see e.g. ). Moreover the need for a more accurate estimate of $`\widehat{g}_Q`$ became recently more evident since the guidance of the chiral extrapolations by the predictions of the heavy meson chiral perturbation theory necessarily involves $`\widehat{g}_Q`$ . In the effective theory based on the combination of the heavy quark and chiral symmetries  one of the essential parameters is the coupling $`\widehat{g}_{\mathrm{}}`$ which appears as the infinite quark mass limit of $`\widehat{g}_Q`$ :

$`\widehat{g}_Q=\widehat{g}_{\mathrm{}}+𝒪\left(1/m_Q^n\right).`$ (3)
While the $`B^{}B\pi `$ decay is forbidden by the lack of phase space and therefore the coupling $`g_{B^{}B\pi }`$ cannot be measured in experiment, the coupling $`g_{D^{}D\pi }`$ has already been measured. CLEO reported recently :

$`g_{D^{}D\pi }`$ $`=`$ $`17.9\pm 0.3\pm 1.9,`$ (4)
$`i.e\widehat{g}_c^{\mathrm{exp}.}`$ $`=`$ $`0.61\pm 0.01\pm 0.07,`$ (5)
where we used eq. (1) to obtain $`\widehat{g}_c`$ from $`g_{D^{}D\pi }`$. Since the pion emerging from this decay is soft, the experimental measurement of $`\mathrm{\Gamma }(D^+)`$ is difficult. The experimental confirmation of the CLEO results would be highly welcome.

Regarding the theoretical estimates the situation is still quite unsatisfactory. Typical values for $`\widehat{g}_Q`$, as obtained by using the QCD sum rules (QCDSR) are low, around $`0.3`$ , in disagreement with $`\widehat{g}_c^{\mathrm{exp}.}`$. Such a disagreement with the experimental value was not, however, observed with similar couplings, such as $`g_{\rho \omega \pi }`$, $`g_{NN\pi }`$, or $`g_{\mathrm{\Sigma }_c^{}\mathrm{\Lambda }_c\pi }`$. To overcome that problem, in ref.  it has been proposed to include the first radial excitations in the hadronic part of the sum rules, below the duality threshold. As a result the coupling $`\widehat{g}_c`$ becomes much larger ($`\widehat{g}_c0.7`$), while the similar modification of the sum rule leaves the decay constants ($`f_D`$ and $`f_D^{}`$) unchanged.

Predictions of various quark models are in the range $`0.3\widehat{g}_{\mathrm{}}0.8`$ . In particular, in the model with the Dirac equation a value of $`\widehat{g}_{\mathrm{}}`$ around 0.6 was found before the experimental result became available . On the lattice, a quenched calculation by UKQCD  led to $`\widehat{g}_{\mathrm{}}=0.42(4)(9)`$. This exploratory study was carried out in the static heavy quark limit where the difficulties to isolate the ground state are well known. In addition it was performed on a coarse lattice, with a low statistics and with only two light quark masses from which it was extrapolated to the chiral limit. Finally, the indirect measurements using the lattice data for the form factors $`F_+`$ and $`F_0`$  have extracted rather small values .

In our recent paper , we presented the results of the quenched lattice simulation in which we used the propagating heavy quarks with the mass around the physical charm quark. We obtained, $`\widehat{g}_c=0.67\pm 0.08_{0.06}^{+0.04}`$, and observed no apparent dependence on the heavy quark mass. Actually the 3 values of $`\widehat{g}_Q`$ measured directly in the neighbourhood of the $`c`$-quark mass were compatible within our error bars and shew only a very small negative slope (see section 4.3 of ref. ). This situation allowed a safe interpolation to the $`c`$-quark but the determination of $`\widehat{g}_{\mathrm{}}`$ relied heavily on the hypothesis that the observed behaviour was valid over the whole mass-range. In order to check the validity of this asumption we have performed a lattice quenched computation of $`\widehat{g}`$ In the static heavy quark limit, the outcome of which is presented in this paper. Our result is

$`\widehat{g}_{\mathrm{}}=0.48\pm 0.03\pm 0.11,`$ (6)
where the errors are statistical and systematic, respectively. Using this and the results we presented in ref. , we were able to interpolate to the $`b`$-quark. We get

$`\widehat{g}_b`$ $`=`$ $`0.58\pm 0.06\pm 0.10,`$ (7)
$`g_{B^{}B\pi }`$ $`=`$ $`47\pm 5\pm 8.`$
By comparing eqs. (6) and (7), our results suggest that the $`1/m_B`$ corrections to $`\widehat{g}_{\mathrm{}}`$ are of the order of $`2030`$%.

The paper is organized as follows: after this introduction we explain in Sec. 1 the relation between $`\widehat{g}`$ and the matrix elements that we compute on the lattice; in Secs. 2 and 3 we present the details of our lattice study and show our main results; Sec. 4 is devoted to the discussion of systematic errors; in Sec. 5 we interpolate to the physical $`B`$-meson mass; finally we summarize our findings in Sec. 6.

## 1 Extracting $`\widehat{g}_{\mathrm{}}`$

The coupling $`g_{H^{}H\pi }`$ (with $`H\{D,B\}`$) is defined by the matrix element

$`H(p^{})\pi (q)|H^{}(p,\lambda )=g_{H^{}H\pi }\left(qϵ^\lambda (p)\right),`$ (8)
where $`q=pp^{}`$ and $`ϵ^\lambda (p)`$ is the polarization of the vector meson. To determine $`g_{H^{}H\pi }`$, we can consider the matrix element of the divergence of the light–light axial current, $`A_\mu =\overline{q}\gamma _\mu \gamma _5q`$. In the limit where the pion is soft, we can write

$`H(p^{})|q_\mu A^\mu |H^{}(p,\lambda )=g_{H^{}H\pi }{\displaystyle \frac{qϵ^\lambda (p)}{m_\pi ^2q^2}}\times f_\pi m_\pi ^2+\mathrm{}`$ (9)
where the dots stand for terms suppressed when $`q^2`$ is small and close to the $`m_\pi ^2`$-pole. To bring out $`g_{H^{}H\pi }`$ from (9) we express the l.h.s in terms of the form factors, $`A_{0,1,2}`$, evaluated on the lattice. The standard parametrization reads

$`H(p^{})\left|A^\mu \right|H^{}(p,\lambda )`$ $`=`$ $`2m_H^{}A_0(q^2){\displaystyle \frac{ϵ^\lambda q}{q^2}}q^\mu +(m_H^{}+m_H)A_1(q^2)\left[ϵ^{\lambda \mu }{\displaystyle \frac{ϵ^\lambda q}{q^2}}q^\mu \right]`$ (10)
$`+A_2(q^2){\displaystyle \frac{ϵ^\lambda q}{m_H+m_H^{}}}\left[p^\mu +p^\mu {\displaystyle \frac{m_H^{}^2m_H^2}{q^2}}q^\mu \right].`$
Note that when $`\stackrel{}{q}=\stackrel{}{0}`$, the soft pion limit is verified, since $`q^2=(m_H^{}m_H)^20`$. Moreover, if we choose $`\stackrel{}{q}=\stackrel{}{0}`$ and $`\stackrel{}{p}=\stackrel{}{p}^{}=\stackrel{}{0}`$, eq. (10) simplifies to

$`H\left|A_i\right|H^{}=(m_H^{}+m_H)A_1(0)ϵ_i^\lambda .`$ (11)
As discussed in our previous paper , when $`q^2=0`$ the form factors $`A_{0,1,2}`$ are related to one another because the $`q^2`$-poles in eq. (10) are unphysical and therefore their residues have to cancel. This leads to the relation

$`2m_H^{}A_0(0)=(m_H^{}+m_H)A_1(0)+(m_H^{}m_H)A_2(0).`$ (12)
When taking the divergence of the axial current in (10), only the term in $`A_0`$ remains on the r.h.s. so that by using eq. (9), we arrive at

$`g_{H^{}H\pi }`$ $`=`$ $`{\displaystyle \frac{2m_H^{}A_0(0)}{f_\pi }}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{f_\pi }}\left[(m_H^{}+m_H)A_1(0)+(m_H^{}m_H)A_2(0)\right].`$ (14)
Finally, inserting the definition (1) leads to the expression for $`\widehat{g}_Q`$

$`\widehat{g}_Q={\displaystyle \frac{(m_H^{}+m_H)}{2\sqrt{m_Hm_H^{}}}}A_1(0)+{\displaystyle \frac{(m_H^{}m_H)}{2\sqrt{m_Hm_H^{}}}}A_2(0).`$ (15)
Owing to the heavy quark symmetry, the vector and the pseudoscalar heavy–light mesons are degenerate in the static limit and eq. (15) simply becomes

$`\widehat{g}_{\mathrm{}}=A_1(0).`$ (16)
Thus, the value of $`\widehat{g}_{\mathrm{}}`$ (in the static heavy quark limit) is simply given by the form factor $`A_1`$ at $`\stackrel{}{q}=\stackrel{}{0}`$.

## 2 Strategy for the calculation of $`\widehat{g}_{_{\mathrm{}}}`$ on the lattice

To determine the form factor $`A_1(0)`$ from the lattice evaluation of the matrix element $`H\left|A^\mu \right|H^{}`$, we use eq. (11) where $`H`$ and $`H^{}`$ now refer to the pseudoscalar and vector mesons containing an infinitely heavy quark. We consider the two- and three-point correlation functions, $`C_2`$ and $`C_{3\mu \nu }`$, with local (L) and smeared (S) interpolating fields for mesons consisting of a static heavy and a propagating light quarks.

In the case in which both the source and the sink are smeared, the two-point functions are defined as

$`C_2^{SS}(t_x)={\displaystyle \underset{\stackrel{}{x}}{}}P^S(x)P_{}^{S}{}_{}{}^{}(0)={\displaystyle \frac{1}{3}}{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \underset{\stackrel{}{x}}{}}V_i^S(x)V_{i}^{S}{}_{}{}^{}(0),`$ (17)
where $`P^S`$ and $`V_\mu ^S`$ are the smeared interpolating fields of the pseudoscalar and vector mesons. $`C_2^{SL}(t_x)`$ is obtained by simply replacing $`P^SP^L=\overline{q}\gamma _5h`$, $`V^SV^L=\overline{q}\gamma _\mu h`$, where $`h`$ stands for the static heavy quark field. The equality between the pseudoscalar and vector meson correlation functions in eq. (17) is a consequence of the heavy quark symmetry. When the excited states are decoupled and the ground state isolated (at large enough $`t_x`$), the two-point function becomes

$`C_2^{SS}(t_x)=(𝒵^S)^2e^{E_H^{SS}t_x},`$ (18)
where $`E_H`$ is the binding energy of the static heavy–light meson, and

$`𝒵^S{\displaystyle \frac{\mathrm{\hspace{0.17em}0}|P^S|H}{(2m_H)^{1/2}}}.`$ (19)
The three-point Green functions with smeared interpolating fields are defined in the following way,

$`C_{3\mu \nu }^{SS}(0,t_x,t_y)={\displaystyle \underset{\stackrel{}{x},\stackrel{}{y}}{}}P^S(y)A_\nu (x)V_{\mu }^{S}{}_{}{}^{}(0)|_{\mathrm{\hspace{0.33em}0}<t_x<t_y},`$ (20)
where the vector and pseudoscalar mesons are inserted at the origin and at a fixed time $`t_y`$, respectively. The insertion time $`t_x`$ of the (local) axial current ($`A_\mu =\overline{q}\gamma _\mu \gamma _5q`$) varies in the range $`0t_y`$. The matrix element we are considering, $`H\left|A_\mu \right|H^{}`$, is then extracted from the following ratio:

$`R^{SS}(t_x)`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{C_{3ii}^{SS}(t_x)𝒵_S𝒵_S}{C_2^{SS}(t_x)C_2^{SS}(t_yt_x)}}`$
$``$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{i=1}{\overset{3}{}}}ϵ_i^\lambda (\stackrel{}{0}){\displaystyle \frac{H\left|A_i(t_x)\right|H^{}}{2m_H}}|_{\mathrm{\hspace{0.33em}0}t_xt_y}=A_1(0)=\widehat{g}_{\mathrm{}}.`$
The time dependence cancels in the ratio, so that a plateau in time of $`R^{SS}(t_x)`$ directly leads to the value of $`A_1(0)`$, which for $`m_q0`$ becomes the desired coupling $`\widehat{g}_{\mathrm{}}`$. In what follows we will use the notation in which $`A_1(0)\widehat{g}_{\mathrm{}}`$, even though this is strictly true only for $`m_q0`$.

In our simulation we choose to work with the usual (improved) Wilson light quarks and with the static heavy quark by using the Eichten–Hill action . In terms of quark propagators and the Wilson line, the three-point correlation function (20) can be written as

$`C_{3\mu \nu }^{LL}(t_x)={\displaystyle \underset{\stackrel{}{x},\stackrel{}{y}}{}}\mathrm{Tr}\left[{\displaystyle \frac{1+\gamma _0}{2}}P_y^0\gamma _\mu S_u(0;x;U)\gamma _\nu \gamma _5S_d(x;y;U)\gamma _5\right]_U,`$ (22)
where $`S_q(x;y;U)=q(x)\overline{q}(y)`$ is the propagator of the light quark $`q`$. In our study the two light quarks are degenerate, $`m_u=m_dm_q`$. The brackets $`\mathrm{}_U`$ indicate the average over the background gauge field configurations ($`U`$), whereas $`P_y^x`$ is the ordered product of links in the time direction (Wilson line),

$`P_y^x=\delta (\stackrel{}{x}\stackrel{}{y}){\displaystyle \underset{t_z=t_x}{\overset{t_y1}{}}}U_t(x+t_z\widehat{t}),`$ (23)
with $`\widehat{t}`$ being the unit vector in the time direction.

To improve the statistical quality of the signal we implement the recent proposal of ref.  and replace the simple link variables in the Wilson line by the so called “fat links”, defined as

$$U_t(x)U_t^{\mathrm{fat}}(x)=\frac{1}{6}\underset{i=x,y,z}{}\left[U_i^{\mathrm{Staple}}(x,x+\widehat{t})+U_i^{\mathrm{Staple}}(x,x+\widehat{t})\right].$$
(24)
In order to isolate the ground state from the correlators (22) at the smallest time separations, it is necessary to devise an efficient smearing procedure. We use the proposal of ref.  and replace the quark fields $`q(x)`$ by

$$q(x)\underset{r=0}{\overset{R_{max}}{}}\left(r+{\scriptscriptstyle \frac{1}{2}}\right)^2\varphi (r)\underset{i=x,y,z}{}\left\{\left[\underset{k=1}{\overset{r}{}}U_i^F(x+\left(k1\right)\widehat{i})\right]q(x+r\widehat{i})+\left[\underset{k=1}{\overset{r}{}}U_i^F^{}(xk\widehat{i})\right]q(xr\widehat{i})\right\}.$$
(25)
The wave function $`\varphi (r)`$ is chosen in such a way that the overlap with the ground state is increased. We take $`\varphi (r)=e^{r/R_b}`$, where $`R_b`$ is a parameter which is fixed by requiring that the smearing be optimal. Note that it is not necessary to normalize the wave function since the normalisation factors cancel in the ratio (2). The smearing also includes the so called fuzzing: the gauge links $`U^F`$ are fuzzed by $`N_F`$ iterations of the following procedure

$$U_j(x)U_j^F(x)=𝒫\left[CU_j(x)+\underset{ij}{}\left[U_i^{\mathrm{Staple}}(x,x+\widehat{j})+U_i^{\mathrm{Staple}}(x,x+\widehat{j})\right]\right],$$
(26)
where $`𝒫`$ is used to project the fuzzed fields onto $`SU(3)`$, and $`C`$ is a real parameter.

## 3 Conditions and results of the simulation

The main results presented in this work are obtained from the simulation performed on a $`24^3\times 28`$ lattice with periodic boundary conditions, at $`\beta =6.2`$. We collect $`160`$ independent $`SU(3)`$ gauge configurations in the quenched approximation. The Wilson fermion action is non-perturbatively $`𝒪(a)`$ improved with $`c_{SW}=1.614`$ . The light quark propagators are computed with the following Wilson hopping parameters:

$`\kappa _q`$ $`=`$ $`0.1344_{q_1},\mathrm{\hspace{0.33em}0.1348}_{q_2},\mathrm{\hspace{0.33em}0.1351}_{q_3},`$ (27)
corresponding to quark masses around the physical strange quark.

We will present our results using our preferred set of parameters, $`t_y=13a`$, $`R_{max}=5a`$, $`R_b=3.0a`$, $`N_F=5`$ and $`C=4`$. The motivation for this choice of parameters, the effect of the smearing (25) and of the modification of the static heavy quark action (24), will be discussed in section 4, where we will also quantify the systematic effects involved in our calculation.

### 3.1 Study of the two-point functions

Light–light operator: The axial current appearing in eq. (20) consists of two light quarks. Even though we use the improved Wilson quarks, for our purpose there is no need to improve the bare axial current since the form factor $`A_1`$ is insensitive to the presence of the $`𝒪(a)`$ improvement term. On the other hand, the improvement of the renormalization constant is of course necessary. We use

$`Z_A^I(g_0^2)=Z_A^{(0)}(g_0^2)\left(1+\stackrel{~}{b}_A(g_0^2)a\rho \right),`$ (28)
where the quark mass $`a\rho `$ is obtained from the axial Ward identity, $`_\mu A_\mu =2\rho P`$ ($`P`$ = $`\overline{q}\gamma _5q`$). The non-perturbatively evaluated renormalization constants are: $`Z_A^{(0)}=0.81(1)`$  and $`\stackrel{~}{b}_A=1.19(6)`$  at $`\beta =6.2`$. We have checked that the inverse lattice spacing obtained from the pion decay constant is consistent with the one we reported in our previous paper. Since the temporal extension of the lattice used in our previous paper is much larger than the present one we will use, whenever needed, $`a^1(f_\pi )=2.71(12)\mathrm{GeV}`$ from ref. .

Static heavy–light correlators: The binding energy $`E_H`$ and the constant $`𝒵`$ are obtained from the two-point correlation function as indicated in eq. (18). To identify the plateau in time where the ground state is properly isolated, we study the effective binding energy,

$`E_{\mathrm{eff}}^{SS}(t)`$ $`=`$ $`\mathrm{ln}\left({\displaystyle \frac{C_2^{SS}(t)}{C_2^{SS}(t+1)}}\right).`$ (29)
The superscript “$`SS`$” refers to the situation in which both interpolating fields are smeared. Note also that the projection $`\frac{1+\gamma _0}{2}`$ of the static quark onto the positive energy states reduces the contamination of the signal coming from the opposite side of the periodic lattice. Therefore, the time reversal symmetry over the time dimension of the lattice is not preserved in the heavy quark limit and the two-point functions are no longer symmetric with respect to $`t=T/2`$. We illustrate in fig. 1 the signals for the effective binding energies of the $`C_2^{SS}`$ and $`C_2^{SL}`$ correlators.

By fitting to the form (18) in the time interval in which the plateau is observed, we extract $`E_H`$ and $`𝒵`$. We have also considered a fit with two exponentials, but did not observe any noticeable effect. The double exponential fit becomes necessary only for when performing simulations without the “fattening” indicated in eq. (24) because in that case the plateaus of $`E_{\mathrm{eff}}(t)`$ are much shorter.

We first calculate the $`B`$-meson decay constant, $`f_{B_s}^{\mathrm{static}}`$. Rather than providing an accurate value for this quantity (a precision determination of $`f_{B_s}^{\mathrm{static}}`$, has already been made in ref. ), our intention is only to test the validity of our computation. $`𝒵^L`$ is not obtained directly from the local–local correlators because the corresponding correlators are much less stable and statistically noisier than the smeared–local ones. For this reason $`𝒵^L`$ is extracted from

$`R_𝒵(t)`$ $`=`$ $`𝒵^S\times {\displaystyle \frac{C_2^{SL}(t)}{C_2^{SS}(t)}}|_{t1}𝒵^L.`$ (30)
In table 1, we present the values of $`E_H^{SS}`$, $`E_H^{SL}`$, $`𝒵^S`$, and $`𝒵^L`$ for a static heavy quark and light quarks with masses corresponding to the $`\kappa _{q_i}`$’s specified in eq. (27).

The decay constant $`f_{B_s}^{\mathrm{static}}`$ is obtained from $`𝒵^L`$ using $`m_{B_s}=5.37`$ GeV, i.e. the physical $`B_s`$-meson mass,

$`f_{B_s}^{static}`$ $`=`$ $`Z_A^{stat}𝒵^L\sqrt{{\displaystyle \frac{2}{m_{B_s}}}}a^{3/2},`$ (31)
and $`Z_A^{stat}=0.77(3)`$, the static heavy–light axial current renormalization constant computed non-perturbatively in ref. . Notice that this value of $`Z_A^{stat}`$ is obtained by using the Eichten–Hill static heavy quark action, i.e. without including the effects of the fattening (24). In ref. , however, it has been argued, on the basis of the observed behaviour of the step scaling function, that the effect on $`Z_A^{stat}`$ of the modification of the static quark action is only very small. Therefore the value we are using is certainly sufficient to achieve our present goal which, we recall, is to check the compatibility of our results with previous ones.

In fig. 2, we provide an illustration of the signal (plateau) of the ratio $`R_z(t)`$, and plot the corresponding $`f_{B_q}^{\mathrm{static}}`$, as obtained by using eq. (31). From the linear interpolation to the strange quark mass (corresponding to $`m_{ss}^2=2m_K^2m_\pi ^2`$), we obtain $`f_{B_s}`$, which we quote as

$`f_{B_s}^{\mathrm{static}}={\displaystyle \frac{Z_A^{stat}}{0.77}}(251\pm 20\mathrm{MeV}).`$ (32)
This value agrees quite well with the one, $`f_{B_s}^{\mathrm{static}}=225\pm 10\mathrm{MeV}`$, recently obtained also at $`\beta =6.2`$, in ref. .

### 3.2 Study of the three-point functions

We observe that, for $`t_x`$ not too close to 0 or $`t_y`$, the three-point function, $`C_{3ii}^{SS}(0,t_x,t_y)`$, introduced in eq. (20), does not depend on $`t_x`$. The reason for this is that, in the static limit, the ground states of the vector and pseudo-scalar mesons are degenerate. Therefore, by writing $`C_{3ii}^{SS}(0,t_x,t_y)`$ as

$`{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}0\left|P(t_y)\right|H_kH_k\left|A_i(t_x)\right|H_l^{}H_l^{}\left|V_{i}^{}{}_{}{}^{}(0)\right|0,`$ (33)
where $`k`$ and $`l`$ label the excitations, one sees that in the diagonal contribution ($`k=l`$) the exponential time dependence writes

$`e^{(E_{H_l^{}}E_{H_k})t_xE_{H_k}t_y}=e^{E_{H_k}t_y},`$ (34)
which is independent of $`t_x`$ ($`t_y`$ is kept fixed). The same is true for the ratio (2) in which the denominator is composed of two correlators $`C_2`$ which are degenerate in this case so that the exponential $`t_x`$ dependence explicitly cancels out, and one ends up with

$`R^{SS}(t_x)={\displaystyle \frac{C_{3ii}^{SS}(t_x)}{𝒵_{S}^{}{}_{}{}^{2}e^{E_H^{SS}t_y}}}.`$ (35)
In other words, the plateau for the ratio $`R^{SS}(t_x)`$ is equivalent to a plateau of $`C_{3ii}^{SS}(t_x)`$. A remarkable feature is that $`C_3`$ already develops a plateau at very a small time $`t_x`$. The physical interpretation could be that in the static limit, the non-diagonal matrix elements of the axial current simply vanish,

$`H_k\left|A_i(t_x)\right|H_l^{}_{kl}0,`$ (36)
which implies that the $`t_x`$-dependence in the l.h.s. of eq. (34) vanishes as well when $`kl`$. Note that this is the case in a quark model framework where the axial current reduces to a light quark spin operator. Since the spatial wave functions of the ground state and the excited states are orthogonal, the spin operator does verify eq. (36222 This statement is verified by our results using the static heavy quark. It is further corroborated by the matrix elements obtained by using the propagating heavy quark with a mass around the physical charm quark: we see that, as the heavy quark mass is increased, the three-point function $`C_{3ii}^{SS}(0,t_x,t_y)`$ begins to develop a plateau for smaller and smaller values of $`t_x`$..

As shown in eq. (2), the value of the bare coupling $`\widehat{g}_{\mathrm{}}^{(0)}`$, is obtained from the fit of $`R^{SS}(t_x)`$ to a constant. We illustrate the signal for $`R^{SS}(t_x)`$ in fig. 3. The corresponding results are listed in table 2, where we also give the final results for $`\widehat{g}_{\mathrm{}}=Z_A^I\widehat{g}_{\mathrm{}}^{(0)}`$, for our set of light quarks (27). We also include the values of $`\widehat{g}_{\mathrm{}}`$ when the light quark is either the strange or the $`u/d`$ quark.

From the linear extrapolation to the chiral limit333We defer to subsection (4.3) the discussion of the validity of this procedure and of the uncertainties it induces., shown in fig. 4, we finally obtain,

$`\widehat{g}_{\mathrm{}}=0.48(3).`$ (37)
The results presented in this section have been obtained after having fixed $`t_y=13a<T/2`$. This choice was made after we performed several simulations with $`t_y/a[8,13]`$ and studied $`R^{SS}(t_x)𝒵_S^2`$ for each $`t_y`$. We have checked that the result for $`R^{SS}(t_x)𝒵_S^2`$ does not change as soon as $`t_y`$ gets larger than $`10a`$. To separate the sources as far as possible we have chosen to fix $`t_y`$ at $`13a`$.

## 4 Smearing, Fattening and Systematic uncertainties

As we already mentioned, to improve the statistical quality of the correlation functions computed with the static heavy quark action, we used the “fattening” procedure . In addition, we use the smearing of the heavy–light interpolating fields , which helps isolating the lowest bound state at lower time separations. In this section we explain how the parameters of the smearing procedure were chosen (25) and show the benefits of using the “fat link” static heavy action instead of the standard Eichten–Hill one. After that discussion, we enumerate the sources of systematic uncertainties and estimate the related errors.

### 4.1 Smearing

The smearing procedure involves four parameters: $`R_{\mathrm{max}}`$, $`R_b`$, $`N_F`$ and $`C`$ (see eqs. (25-26)). The value of $`R_{max}`$ represents half the size of the hypercube where the smearing wave function lives and it is set to the maximal value allowed by the memory limitations. In our case, $`R_{max}=5a`$. The parameter $`R_b`$, which appears in the wave function $`\varphi (r)=e^{r/R_b}`$, is not fixed by any physical requirement other than improving the shape of the plateau for $`R_𝒵`$. The smeared–local correlators are much more sensitive to the variation of $`R_b`$ than the smeared–smeared ones. In fig.5 we show the shape of $`R_𝒵(t)`$ for $`R_b[1.0a,\mathrm{\hspace{0.17em}4.0}a]`$.

We see that if $`R_b3.0a`$, the value of $`R_𝒵(t)`$ on the plateau, and thus also the corresponding $`f_{B_q}^{\mathrm{static}}`$, remains unchanged. The reason is that, as the wave function $`\varphi (r)`$ lives in a typical volume of radius $`R_{max}=5a`$, the shape of the normalized wave function term \[$`(r+1/2)^2\varphi (r)`$\] in eq. (25) does not vary much when changing $`R_b`$ in $`3a5a`$. Of course, all the physical quantities are built in such a way that, at the end, the effect of the normalization of $`\varphi (r)`$ (i.e. the unwanted effect of the smearing) cancels. On the basis of these observations, we fix $`R_b=3a.`$444 As a side remark, we note that the physical Compton wavelength of the $`D`$-mesons $`\lambda _D`$ is given by the mass splitting between the $`D`$ and the $`D_1`$ mesons masses. Its value, $`\lambda _D=2.2\mathrm{fm}`$, is somewhat similar to $`R_b=3.0a2\mathrm{fm}`$, that we use in our simulation.

The remaining two parameters, $`C`$ and $`N_F`$, enter in the fuzzing procedure (26). As far as $`C`$ is concerned, it turns out that its value does not make any significant effect on our result. We varied its value in the range $`C[0.8,\mathrm{\hspace{0.17em}4.0}]`$, and the shape of $`R_𝒵(t)`$ does not change at all (of course, within our present statistical accuracy). On the other hand, the number of fuzzing iterations, $`N_F`$, is chosen to be $`N_F=5`$. The effect of each iteration is to fuzz the spatial links by the gauge fields in the immediate neighborhood. By several iterations we can restore the effect of the wave function inside the volume of radius $`R_b`$, where $`\varphi (r)`$ is supposed to live. As $`N_F`$ iterations cover a volume of typical radius $`(N_F+1)^{1/2}a`$, we should require $`N_F8`$ to fill a volume of radius $`R_b=3.0a`$. It turns out, however, that for $`N_F6`$, the signal for $`R_𝒵`$ begins to increase sharply: a large number of fuzzing iterations destroy the quantum fluctuations as it was observed in the cooling procedures used to study instantons . We checked that $`R_𝒵`$, and thus $`f_B^{\mathrm{static}}`$, remain stable for $`N_F[3,5]`$, and for the final simulation we chose $`N_F=5`$.

### 4.2 Effect of “fattening” the static heavy quark propagator

In this subsection we show that replacing the link variables in the Eichten–Hill action by the “fat” ones (see eq. (24)), does make the plateau for the effective binding energies much larger, including a visible reduction of the statistical noise. In fact, after replacing $`U_tU_t^{\mathrm{fat}}`$ in the Wilson line, the plateau gets extended to $`t12a13a`$, which is actually the upper limit ($`T/2=14a`$). From the comparison of the effective energy, eq. (29), obtained by using the “fat” Wilson line and of the standard one (i.e. without fattening), we can clearly see the improvement in the signal. This is illustrated in fig. 6, where in the case of “fat” Wilson line, the statistics is half what it was in the other case. This improvement is still more striking in the case where one uses the so-called Hyp-fattening (cf. refs. ) instead of the standard one (see ref. ). We see that not only the plateau appears much clearer, but the signal remains longer in time. We also see that the two values of the effective binding energies are different, as expected, since the heavy quark lattice action was modified. In table 3, we compare the results obtained in both cases, i.e. with and without using the fattening procedure.

Before closing this subsection, we note also that, as far as the stability of the signal is concerned, the effect of replacing the standard Eichten–Hill action by the fat-link one is much more rewarding in the case of the heavy–light decay constant than for our $`\widehat{g}_{\mathrm{}}`$. This is so because the operator whose matrix element is related to $`f_B^{\mathrm{static}}`$ includes a static quark, whereas the one related to $`\widehat{g}_{\mathrm{}}`$ does not. The fat link action however helps improving the statistical quality of the signal for $`\widehat{g}_{\mathrm{}}`$. We note that, besides this improvement, the use of fattened links leads to a dramatic decrease of the value of $`f_{B_s}^{\mathrm{static}}`$ with respect to what had been found with the Eichten-Hill action. The fact that the plain EH formalism faces problems when trying to study $`f_B`$ has been known for some time. In refs.  it was shown that the smearing led to systematically lower values555We thank the anonymous referee for drawing those references to our attention.. We show here that a similar effect results from the use of fattened links. In ref.  we show that the Hyp-fattening confirms fully the present results.

### 4.3 Systematic errors

* Smearing:

We combine in quadrature the effects of the variation of all the parameters of the smearing procedure on the resulting value for $`\widehat{g}_{\mathrm{}}`$ (see discussion in sec. 4.1). We estimate that systematics to be $`\pm 15\%`$.
* Discretization errors:

In our study we implemented the full $`𝒪(a)`$ improvement of the Wilson QCD action and the axial current. As we discussed in the text, the improvement of the bare axial current does not influence the value of $`\widehat{g}_{\mathrm{}}`$. As for the renormalization constant, we used the non-perturbatively determined value, including the coefficient $`\stackrel{~}{b}_A`$, which ensures the elimination of the artifacts of $`𝒪(a\rho )`$.

Our main result is obtained from the simulation at $`\beta =6.2`$. In order to study the $`𝒪(a)`$ effects, we also performed the simulations at $`\beta =6.1`$ and $`\beta =6.0`$, keeping the physical volume approximately the same and rescaling the smearing parameters. In fig. 7, we show the chiral behavior of $`\widehat{g}_{\mathrm{}}`$, as computed in all three simulations.

We see that both the results accessed directly from our simulations and those obtained after the linear extrapolation agree within the error bars. To give the reader a more quantitative insight into this effect, we collect in table 4 the results for $`\widehat{g}_{\mathrm{}}`$ and $`f_{B_s}^{\mathrm{static}}`$ from all three simulations.

Since in this work we do not deal with propagating heavy quarks (which are usually the dominant source of discretization errors in the heavy–light meson observables), and since we do not observe any noticeable discretization error on $`\widehat{g}_{\mathrm{}}`$ from our simulations at three different lattice spacings, we will not account for any extra discretization error in our overall systematic uncertainty666Notice also that $`\widehat{g}_{\mathrm{}}`$ is dimensionless so that no discretization errors due to the conversion from the lattice to physical units arise..
* Finite volume:

In this work we did not carry out a detailed study of the finite volume effects. Instead, we will rely on our previous study with the propagating heavy quark where the finite volume effects are estimated to be $`\pm 6`$. This uncertainty will be added in our overall systematic error estimate.
* Chiral extrapolations:

Since we have only three values for the bare light-quark masses, we cannot make a detailed study of the effect of the chiral extrapolation like the one we did in . The values of $`\widehat{g}_{\mathrm{}}`$ that are directly accessed from our simulations follow a very smooth linear behavior when the value of the light quark mass is changed. The result (37) is obtained from the linear extrapolation. To that we will add an error of $`\pm 15\%`$, as estimated in our previous paper . This error is considered symmetric for the reasons explained in detail in ref. : the chiral logs in the full (unquenched) theory have a tendency to lower the value of the coupling $`\widehat{g}_{\mathrm{}}`$, whereas those in the quenched theory do the opposite (increase the value $`\widehat{g}_{\mathrm{}}`$).

After combining the above sources of systematic uncertainty in quadrature, we end up with $`22\%`$ of error to the value given in eq. (37), i.e.

$`\widehat{g}_{\mathrm{}}=0.48\pm 0.03\pm 0.11.`$ (38)
## 5 Heavy mass interpolation to the $`b`$-quark

We can now combine our static heavy quark result, $`\widehat{g}_{\mathrm{}}`$, with those obtained in our previous study in which propagating heavy quarks with masses around the physical charm quark were used , to interpolate to the $`b`$-quark sector. For that purpose we use the spin-averaged mass of the heavy–light mesons, i.e.

$`\overline{m}_H={\displaystyle \frac{3m_V+m_P}{4}}.`$ (39)
Motivated by the heavy quark symmetry, we fit our results to the linear and quadratic forms:

$`\widehat{g}_Q`$ $`=`$ $`\widehat{g}_{\mathrm{}}+{\displaystyle \frac{a_1}{\overline{m}_H}},`$ (40)
$`\widehat{g}_Q`$ $`=`$ $`\widehat{g}_{\mathrm{}}+{\displaystyle \frac{b_1}{\overline{m}_H}}+{\displaystyle \frac{b_2}{\overline{m}_H^2}}.`$ (41)
These fits are shown in fig. 8, from which we read off the values at $`1/\overline{m}_B=(5.314\mathrm{GeV})^1`$, i.e. $`\widehat{g}_b`$:

$`\widehat{g}_b^{(lin)}=0.55(4)(9),\widehat{g}_b^{(quad)}=0.60(7)(9).`$ (42)
As our final estimate we will take the average of the two, and add the difference in the overall systematics. Thus we have

$`\widehat{g}_b=0.58\pm 0.06\pm 0.10,`$ (43)
which by means of eq. (1) gives

$`g_{B^{}B\pi }=47\pm 5\pm 8.`$ (44)
## 6 Summary and conclusions

In this paper we made a study of $`\widehat{g}_{\mathrm{}}`$, the coupling of the lowest-lying doublet of heavy–light mesons to a charged pion, in the infinitely heavy quark limit. The static heavy quark action used in this work is the Eichten and Hill one with fat-links, which ensures a better statistical quality of the correlation functions computed on the lattice. Our result, in the quenched approximation, is

$`\widehat{g}_{\mathrm{}}`$ $`=`$ $`0.48\pm 0.03\pm 0.11,`$ (45)
where the errors are statistical and systematic respectively. The most significant sources of systematic uncertainty are those related to the chiral extrapolation and to the smearing procedures used to suppress the contribution of the higher excited states in the relevant correlation functions.

That result is then combined with the ones reported in our previous work , in which we used propagating heavy quarks with masses close to the one of the physical charm quark. Interpolation to the $`b`$-quark sector leads to

$`\widehat{g}_b=0.58\pm 0.06\pm 0.10\mathrm{and}g_{B^{}B\pi }=47\pm 5\pm 8.`$ (46)
The quantities we have studied in this paper, $`\widehat{g}_{\mathrm{}}`$ and $`g_{B^{}B\pi }`$, are very interesting because

1. they are of high phenomenological relevance,
2. they escape direct measurement,
3. their theoretical determinations are very widespread.

Still, as the attentive reader will have noticed, the errors remain larger large. Referring for instance to eq. (28) above, the statistical error is 6% and the systematic one 22%.

Concerning the statistical error a first factor is obviously the performance level of the available computing facilities. Since this is getting higher and higher rather fast and, at the same time, new computational techniques have been discovered which greatly improve the signal/noise ratio this source of uncertainty will probably be tamed very soon. Some of the systematic errors are no serious problems : most of the various renormalization and improvement constants we need have now been non perturbatively determined with an accuracy at the $`\%`$-level. Note also that some of the (less well-known) heavy-light constants (e.g. $`Z_{B_s}^{static}`$) simply drop out of our computation. Some other ones are more serious :

* The chiral extrapolation (which we have very roughly estimated to contribute by 15%) can be attacked on both the theoretical side (determination of the chiral logs) and the numerical one.
* In this study we have used the static limit only as a limit point for our extrapolations from the finite-mass region. It is possible to go beyond this approach and to determine higher order $`\frac{1}{M_Q}`$ coefficients, which would allow a kind of “educated fit” and to decresase the theoretical uncertainties. This is part of an ambitious program which has been undertaken by the Alpha-group 3 years ago and is still under progress .
* There remains the question of the unquenching. Its possibility while using the fattening procedures has been established (see ref. ) at least in the staggered-fermions case and probably with Wilson fermions too. This is one of the tasks we should undertake in the future.

These are just a few hints regarding what will have to be done now but it is impossible at present to outline any really sensible quantitative perspective regarding the time evolution in the near future of the various uncertainty factors we have just mentionned. Still, pursuing in the direction of a “pure QCD” determination of $`\widehat{g}`$ is clearly of primary importance.

## Acknowledgement

We thank M. Della Morte, J. Reyes and A. Shindler for discussions and comments. The simulation was performed with the APE1000 located in the Centre de Ressources Informatiques (Paris-Sud, Orsay) and purchased thanks to a funding from the Ministère de l’Education Nationale and the CNRS. This work was supported in part by the European Union Human Potential Program under contract HPRN-CT-2000-00/45, Hadrons/Lattice Phenomenology.


format: text

text:
# Dynamic ferromagnetic proximity effect in photoexcited semiconductors

## Abstract

The spin dynamics of photoexcited carriers in semiconductors in contact with a ferromagnet is treated theoretically and compared with time-dependent Faraday rotation experiments. The long time response of the system is found to be governed by the first tens of picoseconds in which the excited plasma interacts strongly with the intrinsic interface between semiconductor and ferromagnet in spite of the existence of a Schottky barrier in equilibrium.

Magnetoelectronics, i.e. the science and technology of using ferromagnets in electronic circuits, is divided into two subfields, metal-based GMR  and semiconductor-based magnetoelectronics Optics , with little common ground.

Metal researchers focus mainly on topics derived from giant magnetoresistance, the large difference in the DC conductance for parallel and antiparallel magnetization configurations of magnetic multilayers, where theoretical understanding has progressed to the stage of materials-specific predictions SSP . More recently, the dynamics of the magnetization vectors in the presence of charge and spin currents has received a lot of attention Slon96 ; Tsoi98 ; Heinrich ; Ralph .

Semiconductor-based magnetoelectronics is motivated by the prospect of integrating new functionalities with conventional semiconductor electronics. The emphasis has been on the basic problem of spin injection into semiconductors, theoretical understanding is less advanced and detailed electronic structure calculations are just starting Wunnicke ; Zwierzycki . Unlike metals, semiconductors can be studied by optical spectroscopies such as the powerful time-resolved Faraday or Kerr rotation techniques, in which a selected component of a spin-polarized excitation cloud in the semiconductor can be monitored on $`\mathrm{ps}`$ time scales Awschalom . In n-doped GaAs, these experiments revealed long spin-coherence of the order of $`\mu \mathrm{s}`$ Awschalom . When the semiconductor is in contact with a ferromagnet, an initially unpolarized electron distribution prepared using linearly polarized light, was found to very quickly acquire a spin polarization - the dynamic ferromagnetic proximity (DFP) effect Kawakami ; Epstein1 ; Epstein2 . In turn DFP efficiently imprints spin information from a ferromagnet onto nuclear spins by dynamic nuclear polarization, opening new options for quantum information storage.

In this Letter we wish to show that metal and semiconductor-based magnetoelectronics can both be understood in terms of coherent spin accumulations. Specifically, the DFP can be treated by the same formalism that successfully describes the dynamics of the magnetization vector in metallic hybrids Gilbert ; Heinrich . The experiments can be understood in terms of a time-dependent interaction between the conduction-band electrons and the ferromagnet in a “fireball-afterglow” scenario. Photoexcited holes are instrumental in helping the electrons to overcome the Schottky barrier between metal and semiconductor in the first $`50\mathrm{ps}`$ (“fireball” regime) and induce the proximity effect. The interaction weakens with vanishing hole density, thus preventing fast decay of the created spin accumulation in the “afterglow”.

Two groups have already contributed important insights into this problem. Ciuti et al. Sham  interpreted the DFP in terms of a spin-dependent reflection of electrons at a ferromagnetic interface through a Schottky barrier in equilibrium, but did not address the time dependence of the problem. Gridnev Gridnev  did investigate the dynamics of the photoexcited carriers, but postulated a phenomenological relaxation tensor with a specific anisotropy that we find difficult to justify. We show how both approaches can be unified and extended by ab initio magnetoelectronic circuit theory Brataas00 ; Bauer03 ; Xia01 ; Xia02 .

We first summarize the experimental evidence Epstein1 ; Epstein2 . Initially, a $`100`$ $`\mathrm{fs}`$ pulse with frequency close to the band gap is absorbed by the semiconductor (100 $`\mathrm{nm}`$ of GaAs). The polarization state is then monitored by time-dependent Faraday rotation measurements of the coherent spin precession in an applied magnetic field. The homogeneously excited carriers (“fireball”) thermalize within a $`\mathrm{ps}`$, in which the holes also lose any initial spin polarization. The interaction time-scale with the ferromagnet (Fe or MnAs) can be deduced from the rise time of the polarization after excitation with a linearly polarized (LP) light pulse to be $`50`$ $`\mathrm{ps}`$ Epstein1 . For long delay times (“afterglow”), the spin relaxation is very slow ($`>2`$ $`\mathrm{ns}`$), comparable to GaAs reference samples in the absence of a ferromagnet. The sample can be also excited by circularly polarized light (CP), in which case the fireball is polarized from the outset. Dynamic nuclear polarization (DNP) by the hyperfine interaction can be detected by deviations of the precession frequency from the bare Larmor frequency, i.e. a modified g-factor Kawakami . DNP should vanish when the external magnetic field is normal to the spin accumulation. It therefore remains to be explained that Epstein et al. Epstein1  observe a modified $`g`$factor for this configuration that differs for LP and CP excitation. Interesting additional information relates to the material dependence, indicating that the polarization induced by Fe is of opposite sign to that induced by MnAs Epstein1 , and to the modulation of the afterglow Larmor frequency by an applied bias Epstein2 .

Let us consider a semiconductor (Sc) film in which a non-equilibrium electron chemical potential $`\mu |=\mu _c,\stackrel{}{\mu }_s|`$ is excited with charge $`\mu _c`$ and spin $`\stackrel{}{\mu }_s|\mu _x,\mu _y,\mu _z|`$ accumulation (in energy units). The bilayer parallel to the $`\left(yz\right)`$ plane consists of a semiconductor in contact with a metallic ferromagnetic film (F) with fixed single-domain magnetization in the direction of the unit vector $`\stackrel{}{m}`$. By its relatively huge density of states a metallic ferromagnet may be treated as a reservoir in equilibrium. A charge $`\left(I_c\right)`$ and spin $`\left(\stackrel{}{I}_s\right)`$ current $`I|=I_c,\stackrel{}{I}_s|`$ (in units of reciprocal time) flows through the ferromagnet/semiconductor (F$`|`$Sc) interface, which is governed by the spin-dependent (dimensionless) conductances $`g_{}`$ and $`g_{}`$ as well as the complex spin-mixing conductance $`g_{}`$ Brataas00 . Physically, the real part of the mixing conductance expresses the angular momentum transfer to and from the ferromagnet, such as the strength of the spin-current induced magnetization torque Slon96 ; Bauer03  or non-local Gilbert damping Gilbert , whereas the imaginary part is an effective magnetic field Huertas00 ; Stiles ; Huertas02 . The microscopic expression for the conductances is Landauer-like

$$g_{\sigma \sigma ^{}}=\underset{nm}{}[\delta _{nm}r_{nm}^\sigma (r_{nm}^\sigma ^{})^{}],$$
(1)
where the reflection coefficient $`r_{nm}^\sigma `$ of an electron in the Sc with spin $`\sigma `$ at the Sc$`|`$F contact between $`n`$th and $`m`$th transverse modes is accessible to ab initio calculations Xia01 ; Xia02 ; Wunnicke ; Zwierzycki . The time dependence of the system is governed by charge and spin conservation Brataas00 :

$`2h𝒟\left({\displaystyle \frac{d\mu _c}{dt}}\right)_{bias}`$ $`=\left(g_{}+g_{}\right)\left(\mu _ce\phi \right)`$
$`+\left(g_{}g_{}\right)\left(\stackrel{}{m}\stackrel{}{\mu }_s\right)`$ (2)
$$2h𝒟\left(\frac{d\stackrel{}{\mu }_s}{dt}\right)_{bias}=2\mathrm{Re}g_{}\stackrel{}{\mu }_s+2\mathrm{Im}g_{}\left(\stackrel{}{m}\times \stackrel{}{\mu }_s\right)$$

$$+\left[\left(g_{}g_{}\right)\left(\mu _ce\phi \right)+\left(g_{}+g_{}2\mathrm{Re}g_{}\right)\left(\stackrel{}{m}\stackrel{}{\mu }_s\right)\right]\stackrel{}{m}$$
(3)
where $`𝒟`$ is the Sc single-spin energy density of states. The electrostatic potential $`\phi `$ due to an applied bias and/or a charge imbalance between electrons and holes will be disregarded in the following (see below). In the presence of a magnetic field $`\stackrel{}{B}`$, the sum of externally applied and hyperfine (Overhauser) fields with ordered nuclear spins, we have to add

$$\left(\frac{d\stackrel{}{\mu }_s}{dt}\right)_{field}=\frac{g_e\mu _B}{\mathrm{}}\stackrel{}{B}\times \stackrel{}{\mu }_s$$
(4)
where $`g_e`$ is the electron g-factor ($`0.4`$ in GaAs) and $`\mu _B`$ the Bohr magneton. These equations can be summarized in terms of a $`4\times 4`$ matrix equation:

$$T_I\frac{d|\mu }{dt}=𝚪|\mu .$$
(5)
where $`T_I=2h𝒟/g`$ is an interface-mediated relaxation time in terms of the total conductance $`g=g_{}+g_{}`$. Choosing $`\stackrel{}{m}`$ parallel to the $`z`$axis:

$$𝚪=\left(\begin{array}{cccc}1& 0& 0& p\\ 0& \eta _r& \eta _i\mathrm{\Omega }_z& \mathrm{\Omega }_y\\ 0& \eta _i+\mathrm{\Omega }_z& \eta _r& \mathrm{\Omega }_x\\ p& \mathrm{\Omega }_y& \mathrm{\Omega }_x& 1\end{array}\right),$$
(6)
$`\left|\mathrm{\Omega }_\alpha \right|/T_I=\left|g_e\right|\mu _BB_\alpha /\mathrm{}`$ is the Larmor frequency, the mixing conductance has been normalized as $`\eta =2g_{}/g`$ with subscripts $`i`$ and $`r`$ denoting its real and imaginary part, and the polarization is defined as $`p=\left(g_{}g_{}\right)/g`$. Eq. $`\left(\text{5}\right)`$ can be solved easily Gridnev  for the boundary conditions corresponding to LP excitation $`\mu ^{LP}\left(0\right)|=1,0,0,0|`$ or circularly polarized (CP) excitation with wave vector in the $`x`$direction $`\mu ^{CP}\left(0\right)|=1,1,0,0|`$.

When the conductance is expressed in terms of an interface transparency parameter $`\kappa `$ times the intrinsic Sc single-spin Sharvin conductance, we can write:

$$T_I\frac{3.5}{\kappa }\frac{L}{100\text{ }\mathrm{nm}}\sqrt{\frac{m^{}}{0.067m_0}\frac{10\text{ }\mathrm{meV}}{\epsilon }}\mathrm{ps}.$$
(7)
For an Ohmic interface with $`\kappa 0.1,`$ $`L_{Sc}=100`$ $`\mathrm{nm}`$ and at a characteristic electron kinetic energy (depending on doping and excitation density) of $`\epsilon =10`$ $`\mathrm{meV},`$ the time constant for GaAs (effective mass $`m^{}=0.067m_0`$) is $`T_I35`$ $`\mathrm{ps}`$, of the order of the experimental rise time of the proximity effect Epstein1 . At long time scales, experiments find $`\stackrel{˘}{T}_I>2`$ $`\mathrm{ns}`$, which corresponds to a strongly reduced transparency of $`\stackrel{˘}{\kappa }0.002.`$

With a few exceptions (notably InAs), Schottky barriers are formed at metal-semiconductor interfaces when interface states in the gap of the semiconductor become filled giving rise to space charges. Photoexcited holes are strongly attracted by the barrier, thereby dragging the electrons with them Flatte  and/or screen the barrier. The observed large $`\kappa `$ in the “fireball” regime reflects the facilitation of electron transport to the Sc$`|`$F interface by the holes. In the limit of predominant ambipolar electron-hole transport (and absence of an applied bias) we can justify negelect of electrostatic potential $`\phi ,`$ but the electron conductance $`g`$ might be affected by the scattering of the holes.

At long time scales, the holes disappear into the ferromagnet or recombination with electrons in the semiconductor, but the electron spin accumulation persists Sham . In this afterglow, remaining space charges vanish when the sample is grounded and net charge transport is suppressed. Eq. (2) thus vanishes and

$$\stackrel{˘}{T}_I\frac{d\stackrel{}{\mu }_s}{dt}=\stackrel{˘}{𝚪}_s\stackrel{}{\mu }_s.$$
(8)
with

$$\stackrel{˘}{𝚪}_s=\left(\begin{array}{ccc}\stackrel{˘}{\eta }_r& \stackrel{˘}{\eta }_i\stackrel{˘}{\mathrm{\Omega }}_z& \stackrel{˘}{\mathrm{\Omega }}_y\\ \stackrel{˘}{\eta }_i+\stackrel{˘}{\mathrm{\Omega }}_z& \stackrel{˘}{\eta }_r& \stackrel{˘}{\mathrm{\Omega }}_x\\ \stackrel{˘}{\mathrm{\Omega }}_y& \stackrel{˘}{\mathrm{\Omega }}_x& 1\stackrel{˘}{p}^2\end{array}\right).$$
(9)
where the decoration $``$ indicates electron transport through the full Schottky barrier.

The kinetic equations are valid when the system is diffuse or chaotic (as a result of interface roughness or bulk disorder). The ferromagnetic elements should have an exchange splitting $`\mathrm{\Delta }`$ which is large enough that the magnetic coherence length $`\mathrm{}_c=\mathrm{}/\sqrt{2m\mathrm{\Delta }}<\mathrm{min}(\mathrm{},L_F),`$ where $`L_F`$ is the thickness of the ferromagnetic layer (typically $`50`$ $`\mathrm{nm}`$). These conditions are usually fulfilled in hybrid systems except for very thin layers with nearly perfect interfaces. The spin relaxation time in GaAs is taken to be very long. We also require that $`L_{Sc}<\sqrt{2\mathrm{}v_F\mathrm{}/3g_e\mu _BB_{ext}}`$ for diffuse systems Huertas00  or $`L_{Sc}<\mathrm{}v_F/g_e\mu _BB_{ext}`$ for ballistic systems, where $`\mathrm{}`$ is the mean free path, $`v_F`$ the Fermi velocity and $`B_{ext}`$ the externally applied field. For samples with $`L_{Sc}=100\mathrm{nm}`$ the applied magnetic fields should therefore not much exceed 5 T.

Thus a microscopic justification can be given for Gridnev’s phenomenological relaxation tensor Gridnev . But whereas we can explain that the longitudinal relaxation $`g`$or $`\stackrel{˘}{g}\left(1\stackrel{˘}{p}^2\right)`$ can differ from the transverse components $`\left(\mathrm{Re}g_{}\text{ or }\mathrm{Re}\stackrel{˘}{g}_{}\right),`$ Gridnev postulated a large difference between the two components *normal* to $`\stackrel{}{m}.`$ Such large magnetic anisotropies can be excluded for Fe and the scenario sketched by Gridnev can not be a generic explanation for all experiments. To explain the proximity polarization, Gridnev’s $`3\times 3`$ Bloch equation *must* be extended to the $`4\times 4`$ kinetic equation (5) that includes a charge current component.

First-principles calculations of the bare interface conductance without Schottky barrier provide a first indication of the transport properties in the first tens of $`\mathrm{ps}`$. We choose here the Fe$`|`$InAs system, which apart from the Schottky barrier, is very similar to Fe$`|`$GaAs, and as such, of great interest in itself. Table 1 summarizes results obtained by scattering matrix calculations with a first-principles tight-binding basis described in Ref. Zwierzycki . We find the reversal of polarization sign with disorder as noted before Zwierzycki , which may explain the negative polarization found for Fe$`|`$GaAs Epstein1 . The real part of the mixing conductance $`\eta _r`$ (torque) is close to unity similar to metallic interfaces, but in contrast to these, the imaginary part $`\eta _i`$ is strongly enhanced. The latter can be explained by the focus on a small number of states with wave vectors close to the origin, which prevents the averaging to zero found in metals Xia02 .

We now model the experiments of Epstein et al. on GaAs$`|`$MnAs (concentrating on Fig. 1 in Epstein1 ) in which an external magnetic field of 0.12 T lies in the interface $`(yz)`$ plane. Analytic solutions of Eq. (5) can be obtained for $`\eta _r=1\left(=\stackrel{˘}{\eta }_r\right)`$(see Table 1 and Xia02 ) such that all modes are exponentially damped by a single interface relaxation time $`T_I.`$ Time is measured in units of $`T_I`$ and the polarization is chosen to be $`p=1\left(=\stackrel{˘}{p}\right)`$. In the fireball regime the quality factor $`\mathrm{\Omega }_y1,`$ and we adopt $`\eta _i=1\left(=\stackrel{˘}{\eta }_i\right)`$. For LP excitation the charge component relaxes in favor of the $`z`$component polarized along the magnetization direction, $`2\mu _z\left(t\right)=e^{\left(1p\right)t/T_I}e^{\left(1+p\right)t/T_I}`$ for $`\mathrm{\Omega }_y=0`$, which is the essence of the DFP effect. The time scale on which the Schottky barrier recovers determines (together with $`p`$) the modulus of the spin accumulation in the afterglow. It is of the same order, but smaller, than $`T_I`$ because of competing electron-hole recombination in the semiconductor. Our Fig. 1 is similar to Fig. 1 in Epstein1 , but additional experimental data on a short time scale are required to guide the development of a more refined model.

As mentioned above, a modified Larmor frequency has been observed Epstein1  even when the photon wave vector is normal to the field. This is at odds with the notion that a CP excited spin accumulation should rotate around the field without net angular momentum transfer to the nuclei. This could be evidence for a DFP effect for the CP configuration. In the fireball, a significant $`\eta _i`$ acts like a magnetic field in the $`z`$direction, causing the initial spin ensemble to precess into the direction of the external magnetic field, which is then able to polarize the nuclear spins. This effect is weaker for LP excitation since, in the brief fireball interval, any spin accumulation has to be generated before it can precess.

For LP excitation the Larmor frequency depends Epstein2  on an applied bias, proving that the electrostatic potential $`\phi `$ can, in general, not be neglected in Eqs. (2,3). This does not invalidate our qualitative arguments since, compared to the bare interface exchange the modifications needed to explain the shifts in the effective Larmor frequencies are small, beyond the accuracy of our model. The decreasing spin lifetime in the afterglow with increasing bias has been explained by inhomogeneus nuclear polarization Epstein2 , but lowering the Schottky barrier by a forward bias also reduces $`\stackrel{˘}{T}_I`$.

In conclusion, we propose a physical picture for the spin dynamics of photoexcited carriers in semiconductor$`|`$ferromagnet bilayers. The experiments can be understood in terms of at least two time scales. In the first 50 ps or so, the photoexcited carriers screen the Schottky barrier efficiently and the interaction of the electrons with the ferromagnet is described by nearly intrinsic interface conductances that can be calculated from first principles. After delay times of $`>`$100 ps, the Schottky barrier protects the semiconductor carriers from fast decay and any residual exchange interaction is very weak. More insight into the interaction of carriers in semiconductors with ferromagnets could be gained by a faster ($`\mathrm{ps}`$) time resolution and higher applied magnetic fields. Quantitative explanation of the experiments requires self-consistent modelling of the combined electron and hole carrier dynamics as well as ab initio calculations of the interface scattering matrices for electrons and holes.

We thank D.D. Awschalom, Y. Kato, and R.J. Epstein for helpful discussions on the Faraday rotation experiments. This work has been supported by the FOM, the NEDO joint research program “Nano Magnetoelectronics”, DARPA award MDA 972-01-1-0024, and the European Commission’s RT Network *Computational Magnetoelectronics* (Contract No. HPRN-CT-2000-00143).


text:
# Dynamical instability of a dark soliton in a quasi-one-dimensional Bose-Einstein condensate perturbed by an optical lattice

## 1 Introduction

The realisation of dilute atomic Bose-Einstein condensates (BEC) in multi-dimensional optical lattices has opened the door to study novel quantum effects within a pure and highly tunable system. These include the quantum phase transition from superfluid state to Mott insulator , interference of coherent matter waves , Bloch oscillations , number squeezed states , and arrays of Josephson junctions . Following the recent experimental creation of dark solitary waves in harmonically trapped BEC  it is a now a realistic prospect to consider these macroscopic excitations in an optical lattice geometry.

A dark soliton is a one-dimensional localised wavepacket consisting of a notch in the ambient density and a phase slip across its centre. It is supported in repulsive/defocussing nonlinear media where the shape of the wavepacket is preserved by a balance between dispersive effects and nonlinearity. Dark solitons have been investigated thoroughly, both experimentally and theoretically, within the analogous field of nonlinear optics . However, the nature of trapped BEC poses new questions for the dynamics and stability of dark solitons. Thermal instabilities, arising from interaction with the thermal cloud, dissipate energy from the soliton, causing it to accelerate to the edge of the condensate and ultimately disappear , while slow solitons are sensitive to quantum fluctuations .

A dark soliton experiences two kinds of dynamical instabilites. Firstly, a dark soliton, which is strictly a one-dimensional object, is prone to transverse instabilities  when embedded in a two or three-dimensional system. For example, dark solitary waves under weak transverse confinement have been observed to bend via the “snake” instability and ultimately decay into vortices, in both BECs  and nonlinear optics . However, by employing tight transverse confinement to freeze out the transverse degrees of freedom, quasi-one-dimensional condensates can be created , in which dark solitons are expected to have greatly enhanced lifetimes . The second dynamical instability arises from the fact that a dark soliton is a solution of a homogeneous system. The longitudinal confinement featured in BEC experiments breaks the integrability of the system, rendering the soliton unstable to motion through the inhomogeneous background density, and leading to decay via sound emission . In previous work we studied this decay mechanism within a harmonic geometry, and isolated two limiting cases. If the emitted sound can escape the system, the power radiated by the dark soliton is proportional to the local soliton acceleration squared , which in turn is intimately linked to the deformation of the soliton profile caused by the emitted sound waves . In the opposite case of an infinite harmonic system, sound reabsorption occurs and the soliton decay becomes fully stabilised. To control this experimentally, we suggested creating a soliton in a tight dimple trap within a weaker harmonic trap. In such a geometry, the depth of the dimple (relative to the chemical potential) provides a sensitive handle on the soliton-sound interactions, which can be observed experimentally by monitoring changes in the soliton amplitude and period .

Recent work by Frantzeskakis et al.  has discussed the dynamics of slowly-moving dark solitons within optical lattices, including the effects of sound emission. The entire range of lattice periodicities were considered in , and, where appropriate, analytic models for the soliton dynamics were derived. In our work we perform detailed quantitative measurements of the sound emission from the soliton and show that the acceleration squared law discussed above is valid (prior to reabsorption of sound). We concentrate on much faster solitons than , and discuss in detail a range of lattice periodicities that can be readily achieved experimentally, where the lattice wavelength is comparable to, or slightly larger than, the soliton width. Given the number of current BEC experiments featuring optical lattices  and the rapid soliton decay that can be induced in certain limits, this system may be preferential to investigate soliton-sound interactions than the dimple trap proposed in .

In Section 2 we introduce the system under consideration, the basic properties of a dark soliton, and the Gross-Pitaevskii equation used to simulate the soliton dynamics. In Section 3 we illustrate the dynamical behaviour of a dark soliton in an optical lattice, in particular the emission of sound, and compare this to the behaviour in the absence of the optical lattice. The dependence of the soliton dynamics on the lattice height and periodicity is discussed in Sections 4 and 5, and in Section 6 we analyse in detail the soliton dynamics and quantify the power radiated. Finally, in Section 7, we make some concluding remarks.

## 2 Theory

The dynamics of a trapped atomic BEC at low temperatures can be well described by the mean-field Gross-Pitaevskii (GP) equation. Assuming a quasi-one-dimensional condensate (oriented along the x-axis), the GP equation reduces to its one-dimensional form,

$`i\mathrm{}{\displaystyle \frac{\psi (x,t)}{t}}=\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{^2}{x^2}}+V(x)+g|\psi (x,t)|^2\mu \right)\psi (x,t).`$ (1)
where $`\psi (x)`$ is the mean-field order parameter, $`m`$ is the atomic mass, and $`\mu =ng`$ is the chemical potential of the effectively one-dimensional system. The nonlinear coefficient $`g`$ arises from the interatomic interactions and has an effective one-dimensional form, $`g=2\mathrm{}^2a/(ml_{}^2)`$, where $`ł_{}=\sqrt{\mathrm{}/(m\omega _{})}`$ is the transverse harmonic oscillator length and $`a`$ is the three-dimensional s-wave scattering length. Repulsive interactions, $`g>0`$, are generally necessary to support a dark soliton, and will be assumed from now on. However, dark soliton-like structures can, under appropriate conditions, be supported by attractive interactions .

We consider an external potential of the form,

$`V(x)={\displaystyle \frac{1}{2}}m\omega _x^2x^2+V_0\mathrm{sin}^2(2\pi x/\lambda ).`$ (2)
The first term represents a harmonic magnetic trap with frequency $`\omega _x`$ and the second term denotes an optical lattice with height $`V_0`$ and periodicity $`\lambda /2`$, where $`\lambda `$ is determined by the wavelength and relative angle between the laser beams forming the lattice. Figure 1 shows a schematic of the trap configuration and initial condensate density.

On a uniform background $`n_0`$, the wavefunction of a dark soliton with speed $`v`$ and position $`(xvt)`$ is given by,

$`\psi (x,t)=\sqrt{n_0}e^{i(\mu /\mathrm{})t}\left(\beta \mathrm{tanh}\left[\beta {\displaystyle \frac{\left(xvt\right)}{\xi }}\right]+i\left({\displaystyle \frac{v}{c}}\right)\right).`$ (3)
where $`\beta =\sqrt{1(v/c)^2}`$ and $`\xi =\mathrm{}/\sqrt{\mu m}`$ is the condensate healing length corresponding roughly to the size of the soliton . The soliton speed $`v`$ depends on its depth $`n_d`$ relative to the background density and the phase slip $`S`$ across its centre via $`v=\sqrt{n_0n_d}=c\mathrm{cos}(S/2)`$, with the maximum speed set by the Bogoliubov speed of sound $`c=\sqrt{\mu /m}`$. The renormalised soliton energy (after subtracting off the background fluid contribution) is given by ,

$`E_s={\displaystyle \frac{4}{3}}\mathrm{}n_0^{3/2}\left[1\left({\displaystyle \frac{v}{c}}\right)^2\right]^{3/2}.`$ (4)
A stationary dark soliton reaches zero density ($`n_d=n_0`$), features an abrupt $`\pi `$ phase slip, and maximum energy, whereas a soliton with maximum speed $`c`$ has no phase or density contrast, and therefore zero energy. Note that, to first order, a dark soliton behaves dynamically like a classical particle .

Throughout this paper we employ harmonic oscillator units, where length is measured in terms of the harmonic oscillator length $`l_x=\sqrt{\mathrm{}/(m\omega _x)}`$ and time in units of $`\omega _x^1`$. In addition a fixed harmonic trap is always considered, corresponding to $`\mu =70\mathrm{}\omega _x`$, and simulations are performed for fixed chemical potential. We will assume realistic experimental parameters: $`\omega _x=2\pi \times 5\mathrm{H}\mathrm{z}`$, $`\omega _{}=250\omega _x`$, and $`\mu _{3D}=8\mathrm{}\overline{\omega }`$, where $`\overline{\omega }=(\omega _x\omega _{}^2)^{1/3}`$. Then, for a $`{}_{}{}^{87}\mathrm{Rb}`$ ($`{}_{}{}^{23}\mathrm{Na}`$) condensate, our space and time units correspond to $`l_x=4.8(9.4)\mu \mathrm{m}`$ and $`\omega _x^1=0.03\mathrm{ms}`$. Furthermore these parameters roughly relate to a healing length $`\xi =0.58(1.13)\mu \mathrm{m}`$, speed of sound $`c=2.5(1.3)\mathrm{mm}/\mathrm{s}`$, one-dimensional peak condensate density $`n_{1D}=2\times 10^7(5\times 10^7)\mathrm{m}^1`$, and total number of atoms $`N=4000(18000)`$. All effects considered in this paper manifest themselves on timescales smaller than the $`1`$s anticipated soliton thermodynamic lifetime .

In order to realistically probe a range of optical lattice periodicities, we envisage employing two almost co-propagating laser beams, of wavelength $`\lambda _{\mathrm{laser}}`$. By varying the angle between the beams $`\theta `$, the effective lattice wavelength can be increased beyond the actual laser wavelength, $`\lambda \lambda _{\mathrm{laser}}/\mathrm{sin}\theta `$. Assuming a laser wavelength $`\lambda _{\mathrm{laser}}0.8\mu \mathrm{m}`$ (typically produced by a diode laser), then lattice wavelengths up to the order of $`\lambda 20\mu \mathrm{m}`$ can easily be produced, corresponding to $`\lambda 5l_x`$ for $`{}_{}{}^{87}\mathrm{Rb}`$.

## 3 Soliton dynamics in an optical lattice and sound emission

We first consider the behaviour of a dark soliton under harmonic confinement only, in order to compare and contrast this to the behaviour observed when an optical lattice is additionally imposed. Figure 2(a) shows a space-time carpet plot of the renormalised condensate density (actual density minus background density) of a dark soliton, initially at the origin with speed $`v=0.5c`$. The soliton oscillates back and forth in the trap at a rate roughly equal to $`\omega _s=\omega _x/\sqrt{2}`$, as expected . As it climbs the side of the trap, passing through the inhomogeneous background density, it emits counter-propagating sound waves of opposing amplitude, and becomes asymmetrically deformed . The sound reflects off the edge of the condensate and subsequently reinteracts with the soliton around $`t5\omega _x^1`$. An equilibrium is rapidly established where the soliton, on average, reabsorbs as much sound as it emits. Hence the soliton energy features oscillations due to the emission/reabsorption process, but with no net decay . The soliton position and energy as a function of time are shown respectively by grey lines in figure 3(a) and (b).

Consider now the case of imposing, in addition to the harmonic confinement, an optical lattice with height $`V_0=0.2\mu `$ and wavelength $`\lambda =2.4l_x`$, as shown in figure 2(b). The presence of the optical lattice causes the soliton dynamics to differ significantly from the pure harmonic case. As the soliton oscillates, it not only climbs the walls of the harmonic trap, but also traverses the lattice bumps. It thereby experiences increased dynamical instability and emits additional sound waves as it passes over each lattice site. The harmonic trap still contributes to the sound emission, as seen by the faint ‘white’ band of sound in figure 2(b), but this emission is heavily obscured (compared to figure 2(a)) by the lattice-generated sound waves. Importantly, in this case, reflection of sound from the edge of the condensate no longer leads to the establishment of an equilibrium between sound emission and reabsortion. This is because the presence of the optical lattice dephases, or mixes, the emitted sound such that only a proportion is reabsorbed. The soliton becomes faster and shallower, and ‘anti-damps’ to larger amplitudes. By $`t20\omega _x^1`$ the soliton has becomes so fast and shallow that it is virtually indistinguishable from, and gradually merges into, the sound field. The soliton position and energy as a function of time are shown respectively by dotted lines in figure 3(a) and (b).

## 4 Effect of varying the lattice height

We now consider the effect on the soliton dynamics of varying the height of the optical lattice $`V_0`$. Throughout this section, the periodicity of the lattice is kept constant at $`\lambda /2=1.2l_x`$. Figure 3 shows the soliton position and energy as a function of time for various lattice heights. The soliton energy $`E_s`$ is calculated by integrating the GP energy functional,

$`\epsilon (\psi )={\displaystyle \frac{\mathrm{}^2}{2m}}|\psi |^2+V|\psi |^2+{\displaystyle \frac{g}{2}}|\psi |^4,`$ (5)
across the ‘soliton region’, conveniently defined to be $`(x_s\pm 5\xi )`$, and subtracting from this the corresponding contribution of the background fluid (obtained from the time-independent solution of equation (1) in the absence of a soliton).

Figure 3 shows soliton oscillations and energies for various lattice heights. Increasing the lattice height leads to two distinct effects which modify the soliton oscillations in opposing manners. Firstly, the optical lattice modifies the condensate profile, periodically reducing the density. This inhibits the motion of the soliton and tends to shift its turning point towards the centre of the trap and increase its oscillation frequency. This effect can be seen by the dashed line ($`V_0=0.1\mu `$) in figure 3(a). Indeed, if the lattice height is sufficiently large, the soliton will initially be trapped within a single lattice site. Secondly, the optical lattice induces sound emission from the soliton. The soliton energy decays essentially monotonically, with the average decay rate increasing with lattice height. This decay leads to a gradual decrease in the oscillation frequency and increase in the amplitude. For small lattice heights, e.g. $`V_0=0.05\mu `$ (long-dashed line in figure 3(a) and (b)) the optical lattice poses little resistance to the soliton motion, and the low level of sound emission it induces causes the main observable effects. These are a small decrease in the soliton oscillation frequency and a slow decay of the soliton energy. The competing effects are most clearly observed for $`V_0=0.4\mu `$ (solid line in figure 3(a)). The soliton is initially confined to a single lattice site, with a high oscillation frequency. Sound emission causes the soliton to gradually ‘anti-damp’ out of this site into adjacent sites, leading to a decrease in its oscillation frequency and an increase in its amplitude. This effect has already been reported in . The soliton rapidly accelerates towards the edge of the cloud where it ultimately merges into the sound field. However, in the limit $`V_0>\mu `$, an array of mini-condensates is formed. Each lattice site becomes an effectively closed system, and the emitted sound cannot escape the site. In this limit the GP equation predicts that full sound reabsorption by the soliton occurs, which stabilises it against decay, as already observed in . In the limit of very small lattice periodicities, the discrete nonlinear Schrodinger equation may be a more appropriate description, in which case additional oscillatory instabilities come into play .

## 5 Effect of varying the lattice periodicity

We next consider how the periodicity of the lattice affects the soliton dynamics, an effect also considered analytically and numerically in . We start by discussing an anomalous regime, the tight lattice regime, where the lattice periodicity is of the order of, or less than, the soliton width. The path and energy of a dark soliton with initial speed $`v=0.5c`$ under the influence of a tight optical lattice ($`\lambda =0.3l_x(2.5\xi )`$ and $`V_0=0.1\mu `$) are shown by the dashed line in figure 4(a) and 4(b) respectively. The path of the soliton is identical to that in the absence of the optical lattice (solid grey line), while the energy has the same form but is shifted to a lower value. This effect can be explained by careful consideration of the density perturbations caused by the optical lattice. In the limit where the lattice wavelength is of the order of, or less than, the soliton width, the fluid cannot heal to the optical lattice, and the associated density modulations smear together. Then the background density becomes reduced across the whole condensate, and in particular the peak density, $`n_0n_0^{}`$. This is due to the fact that, in our simulations, the chemical potential $`\mu `$ is fixed. If, however, the peak density/number of particles were kept fixed, then the chemical potential would increase accordingly. The soliton essentially sees a constant effective optical lattice potential of height $`V_0/2`$, and overall it sees the ambient harmonic potential shifted upwards by this amount, an effect illustrated in figure 4(c) and 4(d). In both cases of a pure harmonic trap and a harmonic trap modified by a tight optical lattice, the initial soliton speed has been set to the same value $`v=0.5c`$ (i.e. the dependence on peak density is removed), ensuring that the paths are identical. If we assume a Thomas-Fermi relationship between potential and density, $`n(x)=n_0V(x)`$, then we would expect, to first order, the peak density to be decreased by $`(V_0/2)n_0=0.05n_0`$. Equation 4 implies that the initial soliton energy will be rescaled by a factor of $`(n_0^{}/n_0)^{\frac{3}{2}}`$, i.e. $`E_s^{\mathrm{init}}0.926E_s^{\mathrm{init}}`$, and this is indeed the modified soliton energy observed.

We now investigate the effect of lattice periodicity on the soliton dynamics. The soliton path and energy for various lattice wavelengths are shown in figure 5(a) and (b) respectively. In order to remove the effect of reduced peak density for sufficiently tight lattices (based on our fixed chemical potential simulations), the soliton energies in figure 5(b) have been rescaled in terms of the initial soliton energy. This effect is quantified in figure 5(c) which shows the variation in the initial soliton energy as a function of lattice wavelength.

Looking now at the decay of the (normalised) soliton energy as a function of time for various lattice wavelengths, we find that increasing the lattice wavelength (but still remaining far from the tight lattice regime considered earlier) leads to the decay of the soliton energy and an associated increase in the soliton amplitude and frequency (e.g. $`\lambda =0.6l_x`$, solid line in figure 5(a) and(b)) in relation to the pure harmonic case (grey line). The rate of this decay increases up to some resonant wavelength ($`\lambda =1.2l_x`$, dashed line). Any further increase in the wavelength leads to a slowing of the decay. This is illustrated in figure 5(b) by the reduced decay rate for $`\lambda =4.8\xi `$ (dotted line) compared to $`\lambda =1.2\xi `$ (dashed line).

The dependence of the soliton decay on the lattice periodicity can be explained in terms of the density inhomogeneity that the soliton experiences. To illustrate this analysis more quantitatively, figure 5(d) plots the loss in soliton energy after one full oscillation as a function of lattice wavelength. For tight lattices, where the lattice periodicity is of the order of, or less than, the soliton width, the soliton effectively sees a displaced harmonic trap, and so experiences no net decay, due to complete sound reabsorption. In this limit, the soliton energy oscillates periodically due to the continuous sound emission/reabsorption process. This takes the form a beating effect caused by the different oscillation frequencies of the soliton and sound (dipole mode) . Therefore, one full oscillation of the soliton in the trap does not correspond to one beat period, and hence the soliton energy has not yet returned to its initial value, explaining the small apparent energy loss. For progressively weaker lattices, the density begins to heal to the optical lattice, thereby modifying the condensate density from the harmonic trap profile. This increases the density inhomogeneity that the soliton sees, and induces an increase in the rate of decay. Furthermore, the presence of these density modulations dephases the sound field such that the soliton is no longer fully stabilised by sound reabsorption.

Increasing the lattice periodicity further leads to some critical wavelength the condensate density can heal fully to the lattice potential. This represents the system with the maximum density inhomogeneity (for a fixed lattice height) and therefore features the most rapid soliton decay. This point is indicated by the resonance in figure 5(d) and occurs for $`\lambda 1.2l_x(10\xi )`$. Further increase in the wavelength does not change the depth of the density modulations, but merely stretches them out. This now reduces the density gradient experienced by the soliton  and causes a reduction in the decay rate. As $`\lambda \mathrm{}`$ the decay will asymptotically tend towards zero as the trap potential returns to becoming effectively harmonic.

## 6 Quantifying the sound emission

We will now quantatively discuss the soliton dynamics up until the time when the sound reinteracts with the soliton. Figure 6 shows the short-term dynamics of a dark soliton in the optical lattice geometry considered in Section 3 ($`V_0=0.2\mu ,\lambda =2.4l_x`$, dotted line of figure 3), along with the corresponding dynamics in the absence of the optical lattice (grey line of figure 3). The soliton motion (figure 6(a)) is clearly inhibited by the modified local density (figure 6(b)), causing it to turn around earlier and at a lower amplitude than in the absence of the optical lattice. The soliton accelerates and decelerates periodically as it traverses the lattice sites, as shown in figure 6(c) (solid line, left axis) with the oscillation weakly modulated by the additional acceleration induced by the ambient harmonic trap (grey line in same figure). This contribution is zero when the soliton is at the centre of the trap and a maximum when the soliton is at an extremum.

The emission of oppositely “charged” sound waves from the soliton leads to an asymmetric distortion of the soliton profile, and a shift of the soliton centre of mass $`x_{\mathrm{cm}}`$ from the density minimum $`x_s`$, where the soliton centre of mass is defined as,

$`x_{\mathrm{cm}}={\displaystyle \frac{_sx(|\psi |^2n)dx}{_s(|\psi |^2n)dx}}.`$ (6)
This shift, shown in figure 6(c) (dashed line, right axis), is directly proportional to the soliton acceleration (solid line, left axis), in accordance with earlier findings .

The motion of the soliton through the lattice causes a periodic cascade of the soliton energy (figure 6(d)). This cascade is more evident when considering the power emitted by the soliton, shown in figure 6(e). The power consists of a series of peaks of various amplitudes, which clearly occur whenever the soliton traverses the side of a lattice site, i.e. where the density gradient is a local maximum. The peaks are greatest at around $`t1.5\omega _x^1`$, which corresponds to the point of maximum soliton acceleration. Indeed, the power emission can be well fitted by an acceleration-squared law,

$`{\displaystyle \frac{dE}{dt}}=\kappa a^2`$ (7)
where $`\kappa `$ is a constant coefficient ($`\kappa =0.0255(\mathrm{}/l_x)`$ in figure 6(e)). This is also in agreement with an acceleration-squared power law derived analytically for dynamically unstable dark solitons in the context of nonlinear optics . Although  deals with a homogeneous system in which the instabilities arise from imperfections in the medium, such an approach has recently been considered for dark solitons in inhomogeneous condensates . When the emitted sound reinteracts with the soliton ($`t>3\omega _x^1`$), it is partially reabsorbed (as evident from the negative values acquired by the power curve of figure 6(e)), and equation (7) ceases to be valid. For the same reasons, the proportionality between acceleration and the deformation parameter is no longer valid.

## 7 Conclusions

We have made a quantitative investigation of the dynamics of a dark soliton induced by the presence of an optical lattice. We find, in general, that the optical lattice significantly affects the motion and dynamical stability of the soliton. The exception to this is in the tight lattice regime where the fluid cannot heal to the lattice and the soliton effectively sees a harmonically trapped condensate.

By considering optical lattices of varying depths and wavelengths we have illustrated the sensitive dependence of the soliton dynamics on these parameters, and explained the dependence in terms of the density inhomogeneity experienced by the soliton. The effect of the optical lattice, and the associated density modulations, is to create extra dynamical instability in the soliton, in addition to that caused by the ambient harmonic trap. This induces soliton decay via the emission of sound waves. For deep and moderately tight lattices this effect is considerable: the soliton rapidly decays, accelerating to the edge of the condensate, where it disappears into the sound field. Unlike in a pure harmonic trap, the presence of the optical lattice tends to modify the emitted sound such that the full stabilisation of the soliton by sound reabsorption is prevented.

The soliton-sound interaction in an optical lattice geometry can be experimentally probed in a quasi-one-dimensional condensate, providing an alternative method to that suggested in . The dynamical dissipation is observable as a change in the soliton amplitude and oscillation frequency, by repeated expansion of the cloud and subsequent absorption imaging, after a variable evolution time in the trap. The timescales considered in this work are consistent with the expected thermodynamic soliton lifetime of 1s . The intensity of the laser beams that form the optical lattice controls the lattice height, while varying the relative angle between the beams allows the full range of lattice periodicities considered here to be realised. Modification of these two parameters enables full experimental control of the dynamical instability of the soliton induced by the optical lattice.

We acknowledge financial support from the UK EPSRC and discussions with D. J. Frantzeskakis and P. G. Kevrekidis.

## References