text:
# Collective Josephson Vortex Dynamics in Long Josephson Junction Stacks

## I Introduction

Dynamics of magnetic vortices in a stack of long Josephson junctions (LJJ) in a magnetic field applied parallel to the junction layers have attracted much attention due to their intriguing applied and fundamental interests. The motion of Josephson vortices in a single junction system has been exploited in various devices. Collective motion of these vortices in a LJJ stack, in which layers of superconductor (S) and insulator (I) are arranged vertically as in Fig. 1, can be exploited in high frequency devices such as tunable submillimeter-wave oscillators and detectors. Here collective motion, including both in-phase and out-of-phase modes shown in Fig. 2, arises from mutual phase-locking of Josephson junctions caused by (magnetic) inductive coupling between screening currents flowing around adjacent Josephson vortices as they move under a bias current. The phase-locking establishes phase coherence across the Josephson junctions. A LJJ stack exhibiting this phase coherence leads to high power output and bandwidth, and it can serve as a model system for scientific studies.

The motion of Josephson vortices in LJJ stacks yields interesting phenomena: (i) Josephson plasma resonance (JPR) and (ii) supercurrent sub-branching. The experiments on both conventional LJJ stacks (e.g., Nb-Al/AlOx-Nb multilayers) and layered superconductors (e.g., Bi2Sr2CaCu2O8+y (BSCCO)) behaving as intrinsic LJJ stacks indicate that JPR can be tuned by magnetic field $`B`$ and can occur over a broad range of frequencies, from microwave to submillimeter-wave. Also, the supercurrent branch in the current-voltage (I-V) data splits into multiple sub-branches when a bias current exceeds some critical value. To explain the data, two theoretical models have been proposed: one is based on the inductive coupling (i.e., magnetic induction model), and the other is based on the coupling due to charge variation in the S layers (i.e., charging effect model).

The magnetic induction model assumes that the S layer thickness $`d_S`$ is much larger than the Debye (charge screening) length $`r_D`$ (i.e., $`d_Sr_D`$), as in conventional LJJ. In this case, charge variations (or electric field) at each S layers are screen out, yielding local charge neutrality. Consequently, the electric field does not lead to the longitudinal coupling between the S layers. In this model, an applied magnetic field induces supercurrents along the S layers and results in the inductive interaction between adjacent S layers. The induction coupling strength is inversely proportional to the common S layer thickness. This model has been used to explain the experimental data for BSCCO. However, the underlying assumption is not justified in BSCCO since $`d_S3\AA `$ and $`r_D23\AA `$.

On the other hand, the charging effect model accounts for the nonequilibrium effect in atomic scale thick superconducting layers. When the S layers are so thin to be comparable to the Debye length (i.e., $`d_Sr_D`$), as in BSCCO, the breakdown of local charge neutrality yields the charging effect. The particle-hole imbalance may also occur since each superconducting layers cannot completely screen out the charge variation. Hence the presence of charge variations yields the interaction between the contiguous superconducting layers and leads to the coupling between the S layers. Recently the charging effect model, neglecting the magnetic induction effect, has been used to interpret the data for BSCCO.

Earlier studies, including numerical simulations of a finite LJJ system, show that these two models can explain the data qualitatively, but considerable inconsistencies between the experimental and the theoretical results have been found. For example, transverse and longitudinal JPR are predicted by the magnetic induction model and the charging effect model, respectively, but the data indicate that both types of resonance occur. Recent experiments on HgI2-intercalated BSCCO and BSCCO single crystals indicate that the supercurrent branch in the I-V data splits into multiple sub-branches in the resistive state when $`BH_o`$ (i.e., low vortex density regime). An estimated value of $`H_o`$ for Nb-Al/AlOx-Nb multilayers and BSCCO is roughly 0.001T and 0.2T, respectively. In the dense vortex regime (i.e., $`BH_o`$), the I-V data exhibit characteristic kinks, and these kinks closely resemble the prediction made by Machida et al., using the magnetic induction model. However a closer examination of the data reveals some inconsistencies. These inconsistencies suggest that a better theoretical model is needed to describe the LJJ stacks.

In this paper, we investigate collective phase dynamics in conventional LJJ stacks and layered superconductors at low magnetic fields (i.e., $`BH_o`$ in which Josephson vortices are in every I layers as in Fig. 2) and at low temperatures (i.e., below the Abrikosov vortex lattice melting temperature), using a theoretical model accounting for both the induction effect and the charging effect. These two effects are equally important in BSCCO since $`r_Dd_S`$, but the charging effect is neglected in many studies because its strength $`\alpha `$ is small (e.g., $`\alpha 0.10.4`$ in BSCCO). We show how the collective motion of Josephson vortices is modified by a weak charging effect. We outline two main results. First, the Josephson plasma dispersion relation, the Swihart velocity, and the stability condition for collective motion in BSCCO are considerably modified from those in Nb-Al/AlOx-Nb multilayers. Second, the splitting of the supercurrent branch in the resistive state is due to collective motion of Josephson vortices, and the width of spread of these sub-branches in the linear I-V regime depends on $`\alpha `$. These results are consistent with the experimental data described above.

The remainder of the paper is organized as follows. In Sec. II, a theoretical model, which accounts for both the magnetic induction effect and the charging effect, is derived by extending previous models. In Sec. III, the Josephson plasma dispersion relation and the Swihart velocity for the collective modes are computed from our model derived in Sec. II. In Sec. IV, we determine the stability condition for the mutually phase-locked modes, performing the linear stability analysis. In Sec. V, we show that the splitting of the supercurrent branch in the resistive state is due to the collective motion of Josephson vortices. Finally, in Sec. VI, we summarize our results and conclude.

## II Theoretical model

In this section, we derive a theoretical model, extending previous approaches. A brief discussion of this model was published. Here we consider a system with a large number of LJJ (i.e., $`N1`$) neglecting the boundary effect and present new results obtained from this model in later sections. To account for both the magnetic induction effect and the charging effect, we start with the gauge-invariant phase difference between the S layers $`\mathrm{}`$ and $`\mathrm{}1`$,

$$\phi _{\mathrm{},\mathrm{}1}=\theta _{\mathrm{}}\theta _\mathrm{}1\frac{2\pi }{\varphi _o}_\mathrm{}1^{\mathrm{}}𝐀𝑑𝐥,$$
(1)
where $`\theta _{\mathrm{}}`$ is the phase of the superconducting order parameter, $`\varphi _o=hc/2e`$ is the flux quantum, and $`𝐀`$ is the vector potential in the I layers. In this paper, we employ the Cartesian coordinates and assume that the S and I layers are stacked along $`z`$-direction and the magnetic field is applied along the $`y`$-direction, as in Fig 1. For simplicity, the thicknesses of the S ($`d_S`$) and I ($`d_I`$) layers are taken to be uniform.

The magnetic induction effect due to the applied magnetic field (along the $`y`$-direction) yields a spatial variation of the phase difference (along the $`x`$-direction). An equation describing the magnetic inductive coupling between the $`S`$ layers

$$\frac{\varphi _o}{2\pi }\frac{\phi _{\mathrm{},\mathrm{}1}}{x}=s(B_{\mathrm{}+1,\mathrm{}}+B_{\mathrm{}1,\mathrm{}2})+d^{}B_{\mathrm{},\mathrm{}1}$$
(2)
is easily obtained by taking a spatial derivative of Eq. (1) and by using the expression for the supercurrent density

$$𝐉_{\mathrm{}}=\frac{\varphi _o}{8\pi ^2\lambda ^2}\left(\theta _{\mathrm{}}\frac{2\pi }{\varphi _o}𝐀_{\mathrm{}}\right).$$
(3)
Here $`d^{}=d_I+2\lambda \mathrm{coth}(d_S/\lambda )`$ and $`s=\lambda [\mathrm{sinh}(d_S/\lambda )]^1`$ are expressed in terms of the London penetration depth $`\lambda `$. The magnetic field $`B_{\mathrm{},\mathrm{}1}`$ in the I layer between two S layers $`\mathrm{}`$ and $`\mathrm{}1`$ is parallel to the layers. Note that $`B_{\mathrm{},\mathrm{}1}`$ differs from $`B`$ since the magnetic field generated by the supercurrent in the S layers modifies the field in the I layer. Using Maxwell’s equation, we express the spatial derivative of the magnetic field as

$$\frac{B_{\mathrm{},\mathrm{}1}}{x}=\frac{4\pi }{c}(J_c\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}J_B+J_{\mathrm{},\mathrm{}1}^T)$$
(4)
where $`J_c`$ is the Josephson critical current density, and $`J_B`$ is a bias current density. Note that the magnetic field entering the I layers yields a triangular Josephson vortex lattice (JVL) when the bias current is either absent or small. The current density $`J^T`$,

$$J_{\mathrm{},\mathrm{}1}^T=\frac{\varphi _o}{2\pi }\frac{\sigma }{𝒟}\frac{\phi _{\mathrm{},\mathrm{}1}}{t}+\frac{ϵ}{4\pi 𝒟}\frac{V_{\mathrm{},\mathrm{}1}}{t},$$
(5)
includes the quasiparticle and the displacement current contribution. Here $`𝒟=d^{}+2s=d_I+2\lambda \mathrm{tanh}(d_S/2\lambda )`$ is the effective thickness of the block layer, $`\sigma `$ is the quasiparticle conductivity, $`ϵ`$ is the dielectric constant of I layers, and $`V_{\mathrm{},\mathrm{}1}`$ is the voltage between the S layers $`\mathrm{}`$ and $`\mathrm{}1`$. Using Eq. (4), we rewrite Eq. (2) as

$`{\displaystyle \frac{\varphi _o}{8\pi ^2}}{\displaystyle \frac{^2\phi _{\mathrm{},\mathrm{}1}}{x^2}}`$ $`=`$ $`J_c\left[s(\mathrm{sin}\phi _{\mathrm{}+1,\mathrm{}}+\mathrm{sin}\phi _{\mathrm{}1,\mathrm{}2})+d^{}\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}\right]𝒟J_B`$ (6)
$`+`$ $`d^{}J_{\mathrm{},\mathrm{}1}^T+s(J_{\mathrm{}+1,\mathrm{}}^T+J_{\mathrm{}1,\mathrm{}2}^T).`$ (7)
Note that $`𝒮=s/d^{}`$ measures the induction coupling strength, and $`𝒮=0.5`$ in the strong coupling limit. The phase difference equation (Eq. (4) in Ref. 27) derived by Bulaevskii and Clem within the framework of Lawrence-Doniach model can be obtained from Eq. (7) when the time-dependent terms are neglected (i.e., $`J_{\mathrm{},\mathrm{}1}^T=0`$) and the relations $`(8\pi ^2/\varphi _o)J_cd^{}=(2/\lambda _J^2)+(1/\lambda _c^2)`$ and $`(8\pi ^2/\varphi _o)J_cs=1/\lambda _J^2`$ are used. Here $`\lambda _c`$ is the magnetic penetration depth in the direction perpendicular to the S layers.

The presence of a nonequilibrium state leads to the interaction between the S layers. When the S layer thickness is comparable to the Debye screening length (i.e., $`r_Dd_S`$), the S layers are in a nonequilibrium state because the charge variations in these layers are not completely screened. This incomplete charge screening enhances the temporal variation of the phase difference. One can include this effect in the phase dynamics, modifying the usual AC Josephson relation, which is a time derivative of Eq. (1), to

$$\frac{\varphi _o}{2\pi }\frac{\phi _{\mathrm{},\mathrm{}1}}{t}=V_{\mathrm{},\mathrm{}1}+\mathrm{\Phi }_{\mathrm{}}\mathrm{\Phi }_\mathrm{}1$$
(8)
as a way to account for a nonzero gauge-invariant potential $`\mathrm{\Phi }_{\mathrm{}}=\varphi _{\mathrm{}}+(\mathrm{}/2e)(\theta _{\mathrm{}}/t)`$ generated inside the S layers. Here $`\varphi _{\mathrm{}}`$ is the electrostatic potential. The modified AC Josephson relation of Eq. (8) can be rewritten as

$$\frac{\varphi _o}{2\pi }\frac{\phi _{\mathrm{},\mathrm{}1}}{t}=V_{\mathrm{},\mathrm{}1}\alpha (V_{\mathrm{}+1,\mathrm{}}2V_{\mathrm{},\mathrm{}1}+V_{\mathrm{}1,\mathrm{}2})\mathrm{\Psi }_{\mathrm{}}+\mathrm{\Psi }_\mathrm{}1,$$
(9)
using the charge density $`\rho _{\mathrm{}}=(\mathrm{\Phi }_{\mathrm{}}\mathrm{\Psi }_{\mathrm{}})/4\pi r_D^2`$ and the Maxwell’s equation $`ϵ𝐄=4\pi \rho `$. Here $`\alpha =ϵr_D^2/𝒟d_S`$ measures the strength of the charging effect, and $`\mathrm{\Psi }_{\mathrm{}}`$ measures the particle-hole imbalance in the S layer. For simplicity, we consider only the charging effect by setting $`\mathrm{\Psi }_{\mathrm{}}=\mathrm{\Psi }_\mathrm{}1=0`$, as it has been done in earlier studies. Note that the usual AC Josephson relation is obtained from Eq. (9) when $`\alpha =0`$. This indicates that the charging effect (i.e., $`\alpha 0`$) enhances the coupling between neighboring junctions. Using Eq. (5), we relate the time derivative of the phase difference to the current densities and obtain

$`J_c\left[{\displaystyle \frac{1}{\omega _p^2}}{\displaystyle \frac{^2\phi _{\mathrm{},\mathrm{}1}}{t^2}}+{\displaystyle \frac{\beta }{\omega _p}}{\displaystyle \frac{\phi _{\mathrm{},\mathrm{}1}}{t}}{\displaystyle \frac{\alpha \beta }{\omega _p}}{\displaystyle \frac{}{t}}\left(\phi _{\mathrm{}+1,\mathrm{}}2\phi _{\mathrm{},\mathrm{}1}+\phi _{\mathrm{}1,\mathrm{}2}\right)\right]=`$ (10)
$`J_{\mathrm{},\mathrm{}1}^T\alpha (J_{\mathrm{}+1,\mathrm{}}^T2J_{\mathrm{},\mathrm{}1}^T+J_{\mathrm{}1,\mathrm{}2}^T).`$ (11)
Here $`\omega _p=c/(\sqrt{ϵ}\lambda _c)`$ is the plasma frequency, $`\beta =(4\pi /c)(\sigma \lambda _c/\sqrt{ϵ})=1/\sqrt{\beta _c}`$, and $`\beta _c`$ is the McCumber parameter. Note that the spatial variation of $`\phi _{\mathrm{},\mathrm{}1}`$ is neglected in the charging effect model of Eq. (11). The terms of the order $`𝒪(\alpha \beta )`$ can be safely neglected since $`\alpha `$ and $`\beta `$ are small in the layered superconductors (i.e., $`\alpha \beta 1`$). For example, the experimental value for $`\alpha `$ and $`\beta `$ in BSCCO are roughly 0.1-0.4 and 0.2, respectively. Neglecting these small terms, we rewrite Eq. (11) as

$$J_c\left(\frac{1}{\omega _p^2}\frac{^2\phi _{\mathrm{},\mathrm{}1}}{t^2}+\frac{\beta }{\omega _p}\frac{\phi _{\mathrm{},\mathrm{}1}}{t}\right)J_{\mathrm{},\mathrm{}1}^T\alpha (J_{\mathrm{}+1,\mathrm{}}^T2J_{\mathrm{},\mathrm{}1}^T+J_{\mathrm{}1,\mathrm{}2}^T).$$
(12)
As we shall see in Sec. III, the charging effect terms in Eq. (12) yield purely longitudinal Josephson plasma excitations.

A theoretical model, including both the magnetic induction effect and the charging effect, can be obtained easily by noting that the magnetic induction model of Eq. (7) and the charging effect model of Eq. (12) are coupled to each other via the current density $`J^T`$. Combining Eqs. (7) and (12), we obtain the coupled sine-Gordon equations,

$`{\displaystyle \frac{^2\phi _{\mathrm{},\mathrm{}1}}{x^2}}{\displaystyle \frac{1}{\lambda _c^2𝒟}}\left({\displaystyle \frac{1}{\omega _p^2}}{\displaystyle \frac{^2\mathrm{\Delta }_{\mathrm{}}}{t^2}}+{\displaystyle \frac{\beta }{\omega _p}}{\displaystyle \frac{\mathrm{\Delta }_{\mathrm{}}}{t}}\right){\displaystyle \frac{\alpha }{\lambda _c^2𝒟}}{\displaystyle \frac{1}{\omega _p^2}}{\displaystyle \frac{^2\mathrm{\Xi }_{\mathrm{}}}{t^2}}=`$ (13)
$`{\displaystyle \frac{1}{\lambda _c^2𝒟}}\left[d^{}\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}+s(\mathrm{sin}\phi _{\mathrm{}+1,\mathrm{}}+\mathrm{sin}\phi _{\mathrm{}1,\mathrm{}2})𝒟{\displaystyle \frac{J_B}{J_c}}\right],`$ (14)
where $`\mathrm{\Delta }_{\mathrm{}}=d^{}\phi _{\mathrm{},\mathrm{}1}+s(\phi _{\mathrm{}+1,\mathrm{}}+\phi _{\mathrm{}1,\mathrm{}2})`$ and $`\mathrm{\Xi }_{\mathrm{}}=s(\phi _{\mathrm{}2,\mathrm{}3}+\phi _{\mathrm{}+2,\mathrm{}+1})+(d^{}2s)(\phi _{\mathrm{}1,\mathrm{}2}+\phi _{\mathrm{}+1,\mathrm{}})+2(sd^{})\phi _{\mathrm{},\mathrm{}1}`$. The third term on the left hand side of Eq. (14) (due to the charging effect) is the main modification from the earlier models. Hence, Eq. (14) becomes identical to the phase difference equation derived by Bulaevskii et al. (Eq. (11) in Ref. 17) when $`\alpha =0`$. Using Eq. (14), we show below that a weak charging effect in BSCCO (i.e., $`\alpha 0.10.4`$) can yield significant changes to the phase dynamics.

## III Josephson plasma dispersion relation

We now determine the dispersion relation for the Josephson plasma and the Swihart velocity for the collective modes, using linear analysis: $`\phi _{\mathrm{},\mathrm{}1}=\phi _{\mathrm{},\mathrm{}1}^{(0)}+\phi _{\mathrm{},\mathrm{}1}^{}`$. Here $`\phi _{\mathrm{},\mathrm{}1}^{}`$ describes small fluctuations about $`\phi _{\mathrm{},\mathrm{}1}^{(0)}`$ describing uniform motion of Josephson vortices in the I layer between $`\mathrm{}`$ th and $`\mathrm{}1`$ th S layers. $`\phi _{\mathrm{},\mathrm{}1}^{(0)}`$ is zero in the Meissner state, but in general, it depends on a magnetic field, allowing JPR to be tunned by the field. The effect of magnetic field on JPR can be accounted for more accurately via the field dependence of $`J_c`$ and via imposing the boundary condition, $`(\varphi _o/2\pi )(\phi _{\mathrm{},\mathrm{}1}/x)=𝒟B`$, explicitly at both $`x=0`$ and $`x=L_x`$ (a junction length) in numerical simulations.

When the bias current $`J_B`$ equals the Josephson current in each of the I layers (i.e., $`J_B=J_c\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}^{(0)}`$), we describe the motion of vortices in terms of a uniform motion $`\phi _{\mathrm{},\mathrm{}1}^{(0)}`$ and small perturbation $`\phi _{\mathrm{},\mathrm{}1}^{}`$ about $`\phi _{\mathrm{},\mathrm{}1}^{(0)}`$. The uniform phase motion is described by

$$\frac{^2\phi _{\mathrm{},\mathrm{}1}^{(0)}}{x^2}\frac{1}{\lambda _c^2𝒟}\left(\frac{1}{\omega _p^2}\frac{^2\mathrm{\Delta }_{\mathrm{}}^{(0)}}{t^2}+\frac{\beta }{\omega _p}\frac{\mathrm{\Delta }_{\mathrm{}}^{(0)}}{t}\right)\frac{\alpha }{\lambda _c^2𝒟}\frac{1}{\omega _p^2}\frac{^2\mathrm{\Xi }_{\mathrm{}}^{(0)}}{t^2}=0,$$
(15)
while small fluctuations (i.e., $`\phi _{\mathrm{},\mathrm{}1}^{}`$) about $`\phi _{\mathrm{},\mathrm{}1}^{(0)}`$ are described by

$`{\displaystyle \frac{\varphi _o}{8\pi ^2}}{\displaystyle \frac{^2\phi _{\mathrm{},\mathrm{}1}^{}}{x^2}}=J_c[d^{}\phi _{\mathrm{},\mathrm{}1}^{}\mathrm{cos}\phi _{\mathrm{},\mathrm{}1}^{(0)}+s(\phi _{\mathrm{}+1,\mathrm{}}^{}\mathrm{cos}\phi _{\mathrm{}+1,\mathrm{}}^{(0)}+\phi _{\mathrm{}1,\mathrm{}2}^{}\mathrm{cos}\phi _{\mathrm{}1,\mathrm{}2}^{(0)})]`$ (16)
$`+d^{}J_{\mathrm{},\mathrm{}1}^T+s(J_{\mathrm{}+1,\mathrm{}}^T+J_{\mathrm{}1,\mathrm{}2}^T),`$ (17)
$`J_c\left({\displaystyle \frac{1}{\omega _p^2}}{\displaystyle \frac{^2\phi _{\mathrm{},\mathrm{}1}^{}}{t^2}}+{\displaystyle \frac{\beta }{\omega _p}}{\displaystyle \frac{\phi _{\mathrm{},\mathrm{}1}^{}}{t}}\right)=J_{\mathrm{},\mathrm{}1}^T\alpha (J_{\mathrm{}+1,\mathrm{}}^T2J_{\mathrm{},\mathrm{}1}^T+J_{\mathrm{}1,\mathrm{}2}^T).`$ (18)
Equations (17) and (18) are coupled through the current density $`J_{\mathrm{},\mathrm{}1}^T`$, suggesting that these equations can be simplified by expressing $`\phi _{\mathrm{},\mathrm{}1}^{}`$ and $`J_{\mathrm{},\mathrm{}1}^T`$ as Fourier series in space for the $`z`$-direction: $`\phi _{\mathrm{},\mathrm{}1}^{}=_{m=1}^{2N+1}T_me^{ik_m\mathrm{}a}`$ and $`J_{\mathrm{},\mathrm{}1}^T=_{m=1}^{2N+1}J_m^Te^{ik_m\mathrm{}a}`$. $`k_m=m\pi /(N+1)a`$ represents the wavenumber for the collective mode along the $`z`$-direction, $`a=d_I+d_S`$, $`m`$ is the mode index, and $`N`$ represents the number of Josephson junctions in a stack.

The Josephson plasma mode dispersion relation is determined easily by approximating that $`\phi _{\mathrm{},\mathrm{}1}^{(0)}\phi ^{(0)}`$ and by combining Eqs. (17) and (18) into a single equation as

$$\frac{^2T_m}{t^2}+\omega _p^2\left[(k_x\lambda _c)^2𝒜_m_m+𝒜_m\mathrm{cos}\phi ^{(0)}\right]T_m0,$$
(19)
where $`𝒜_m=1+4\alpha \mathrm{sin}^2(k_ma/2)`$, $`_m=\{1+4[𝒮/(1+2𝒮)]\mathrm{sin}^2(k_ma/2)\}^1`$, and $`𝒮=s/d^{}`$. Here, we set $`\beta =0`$ for simplicity. From Eq. (19), we obtain the dispersion relation of

$$\omega (k_x,k_m)=\omega _p\left[(k_x\lambda _c)^2𝒜_m_m+𝒜_m\mathrm{cos}\phi ^{(0)}_t\right]^{1/2}$$
(20)
for the collective mode. $`\mathrm{}_t`$ represents thermal averages. The dispersion relation of Eq. (20) naturally recovers both purely longitudinal and purely transverse plasma excitations at $`k_x=0`$ and at $`k_m=0`$, respectively. However, there are notable differences between our result of Eq. (20) and the results from other models. Figure 3 illustrates the difference between the dispersion relation of our model and that of the magnetic induction model (Fig. 3(a)) and that of the charging effect model (Fig. 3(b)).

The changes in the dispersion relation due to the charging effect increase the characteristic velocity of the collective mode. The group velocity for the electromagnetic waves in these LJJ is easily determined from Eq. (20) by evaluating

$$\frac{d\omega }{dk_x}=\frac{\omega _p^2}{\omega }k_x\lambda _c^2𝒜_m_m$$
(21)
within the linearized model. This group velocity, asymptotically (i.e., as $`k_x\mathrm{}`$), leads to the Swihart velocity

$$\overline{c}_m=c_o\left[\frac{1+4\alpha \mathrm{sin}^2(k_ma/2)}{1+4\left(\frac{𝒮}{1+2𝒮}\right)\mathrm{sin}^2(k_ma/2)}\right]^{1/2},$$
(22)
the effective maximum velocity for the collective mode $`m`$. Here $`c_o=c/\sqrt{ϵ}`$. Equation (22) recovers the result of the magnetic induction model ($`\overline{c}_m^{MI}`$) when $`\alpha =0`$ (i.e., $`\overline{c}_m^{MI}=\overline{c}_m(\alpha =0)`$), indicating that the charging effect yields the mode-dependent enhancement of the Swihart velocity from $`\overline{c}_m^{MI}`$. For example, the Swihart velocity is not enhanced for the $`m=1`$ mode (i.e., $`\overline{c}_1=\overline{c}_1^{MI}`$), but it is enhanced for the $`m=N`$ mode (i.e., $`\overline{c}_N=(1+4\alpha )^{1/2}\overline{c}_N^{MI}`$). This enhancement reflects the increase in the coupling strength between the junctions due to the charging effect and indicates that the threshold velocity $`v_{th}`$ (=$`\overline{c}_N`$) for emitting Cherenkov radiation (i.e., non-Josephson emission) is also increased. For example, $`v_{th}=\overline{c}_N^{MI}`$ when $`\alpha =0`$, but $`v_{th}=1.34\overline{c}_N^{MI}`$ when $`\alpha =0.2`$. Evidence, indicating the need to account for the charging effect, may be also found in the I-V data for BSCCO. Recent analysis of the I-V data in the dense vortex regime (i.e., $`BH_o`$) indicates that a better agreement between the predicted and observed position of the kinks can be obtained if the Swihart velocity for $`m>1`$ is slightly larger than $`\overline{c}_m^{MI}`$. This suggests that accounting for the charging effect is important for quantitative understanding of the kinks in the I-V curves.

In Fig. 4, we compare the Josephson plasma mode dispersion for (a) Nb-Al/AlOx-Nb multilayers and (b) BSCCO in the Meissner state (i.e., $`\phi ^{(0)}_t0`$), using of Eq. (20). To illustrate the difference between the dispersion of collective mode for these two systems, we use the experimental values for the parameters $`\alpha `$ (i.e., charging effect strength) and $`𝒮`$ (i.e., induction coupling strength). For the spectrum corresponding to the Nb-Al/AlOx-Nb multilayers (Fig. 4(a)), we chose $`\alpha =0.0`$ and $`𝒮0.47`$ (assuming $`\lambda 900\AA `$, $`d_I20\AA `$, and $`d_S30\AA `$). Here we chose $`\alpha =0`$ since the charging effect is negligible when $`d_S`$ is much larger than an atomic length. For the spectrum corresponding to BSCCO (Fig. 4(b)), we chose $`\alpha =0.2`$ and $`𝒮0.49999`$ (since $`\lambda 1500\AA `$, $`d_I15\AA `$, and $`d_S3\AA `$). Here $`\alpha =0.2`$ is chosen. There are two notable differences between Figs. 4(a) and 4(b). First, due to a stronger inductive coupling (i.e., $`𝒮=0.49999`$ versus $`0.47`$), the frequency $`\omega /\omega _p`$, for a fixed $`k_x\lambda _c`$, near $`k_ma=0`$ in Fig. 4(b), decreases more sharply with $`k_ma`$ than that in Fig. 4(a). Second, due to the charging effect (i.e., $`\alpha =0.2`$ versus 0.0), the collective mode frequency for $`k_x\lambda _c=0`$ shows a dispersion as a function of $`k_m`$ in Fig. 4(b), indicating purely longitudinal excitations, while no dispersion is shown in Fig. 4(a), indicating the absence of these excitations. Note that the effect of finite, but small, $`\beta `$ is negligible, here.

## IV Stability of collective modes

In this section, we discuss the stability of uniform motion of collective modes (i.e., moving JVL) shown in Fig. 2 against small fluctuations. The structure of the moving JVLs, driven by a bias current, evolves as a function of its velocities. This evolution can be easily understood in terms of the stable-unstable transition for the collective modes.

We now carry out the linear analysis and determine the condition for maintaining stable uniform motion (i.e., the condition for bound oscillations of $`\phi _{\mathrm{},\mathrm{}1}^{}`$) by computing the velocities at which the driven collective modes are stable. Here, instability of uniform motion arises when the amplitude of fluctuations grows exponentially as the collective modes propagate along the junction layers. Similar analysis, not including the charging effect, have been carried out to investigate the stability of moving JVL against lattice deformation. Also the effects of quantum and thermal fluctuations have been studied. We note that accounting for either the dynamic phase transition induced by the lattice displacements or the fluctuation effects are beyond the scope of the present analysis.

We proceed the analysis writing the spatial and the temporal dependence of phase fluctuations (i.e., $`\phi _{\mathrm{},\mathrm{}1}^{}`$) of Eqs. (17) and (18) in Fourier space for the $`z`$-direction. Combining Eqs. (17) and (18) in Fourier space, we obtain

$$\frac{\lambda _c^2𝒟}{d^{}𝒞_m}\frac{^2T_m}{x^2}\frac{𝒜_m}{\omega _p^2}\frac{^2T_m}{t^2}\left(\frac{\beta }{\omega _p}\frac{T_m}{t}+\mathrm{cos}\phi _m^{(0)}T_m\right)=0$$
(23)
where $`𝒞_m=1+2𝒮\mathrm{cos}k_ma`$. The uniform motion of the phase locked mode $`\phi _m^{(0)}`$ with the wavenumber $`k_m`$ is given by $`\phi _m^{(0)}=\overline{\omega }_mt+k_xx+\phi _{o,m}`$, where $`\overline{\omega }_m=(2\pi /\varphi _o)V_m𝒜_m`$ is the Josephson frequency, $`V_m`$ is the average voltage, and $`\phi _{o,m}`$ is a mode dependent constant. Here, the induced field contribution to $`\phi _m^{(0)}`$ from the Josephson effect is neglected. This contribution is negligible when the magnetic vortices are in every I layers, as the case for $`BH_o`$. $`\overline{\omega }_m/k_x`$ is the velocity of the collective mode $`m`$. We transform Eq. (23) into a familiar Mathieu equation in the following two steps: first, make a change of variables from $`(x,t)`$ to $`\zeta _m(=\phi _m^{(0)})`$; and second, let $`T_m=\overline{T}_me^{\mathrm{\Gamma }_m\zeta _m/2}`$. Here $`\mathrm{\Gamma }_m=\beta \overline{\omega }_m\omega _p𝒞_m/\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_m=𝒜_m𝒞_m\omega _m^2(k_x\lambda _c)^2(𝒟/d^{})\omega _p^2`$. The stability condition is determined, solving

$$\frac{^2\overline{T}_m}{\zeta _m^2}+[\overline{\delta }_m^T+\overline{\eta }_m^T\mathrm{cos}\zeta _m]\overline{T}_m=0,$$
(24)
where $`\overline{\delta }_m^T=\mathrm{\Gamma }_m^2/4`$ and $`\overline{\eta }_m^T=𝒞_m\omega _p^2/\mathrm{\Omega }_m`$ are the mode dependent (i.e., $`k_m`$) parametric constants. Solutions of Eq. (24) exhibit instability for certain values of $`\overline{\delta }_m^T`$ and $`\overline{\eta }_m^T`$, indicating that the collective mode becomes unstable against small fluctuations. We determine the stability condition, finding $`(\overline{\delta }_m^T,\overline{\eta }_m^T)`$ at which all solutions of Eq. (24) are bounded. Note that a similar parametric instability (in the $`\overline{\delta }_m^T\overline{\eta }_m^T`$ space) occurs both in the magnetic induction model (i.e., $`\alpha =0`$) and in the charging effect model (i.e., $`k_x\lambda _c=0`$).

For finding the stability condition for the collective modes, it is useful to determine, first, the stability diagram of the Mathieu equation

$$\frac{d^2T}{d\zeta ^2}+[\delta +\eta \mathrm{cos}\zeta ]T=0,$$
(25)
and then, find the values of $`(\overline{\delta }_m^T,\overline{\eta }_m^T)`$ corresponding to the stable region of this diagram. The boundary curves separating the region of bound (stable) and unbound (unstable) solutions can be obtained easily, solving Eq. (25) numerically following the procedure outlined in Ref. 35. The boundary curves for the periodic oscillations with the period $`2\pi `$ (i.e., $`\zeta =2\pi `$) and $`4\pi `$ (i.e., $`\zeta =4\pi `$) are obtained, imposing that the determinant $`_n`$ for $`n=\mathrm{}`$, derived from Eq. (25) writing $`T=_{n=\mathrm{}}^{n=\mathrm{}}C_ne^{in\zeta }`$ for $`\zeta =2\pi `$ and $`T=_{n=\mathrm{}}^{n=\mathrm{}}d_ne^{in\zeta /2}`$ for $`\zeta =4\pi `$, is zero (i.e., $`_n=0`$). Here $`_n`$ is the determinant of a $`(2n+1)\times (2n+1)`$ matrix for a periodic solution with the period $`2\pi `$ (or a $`2n\times 2n`$ matrix for a periodic solution with the period $`4\pi `$). The determinant $`_n`$ can be computed using the recursion relation

$$_{n+2}=(1\gamma _{n+2}\gamma _{n+1})_{n+1}\gamma _{n+2}\gamma _{n+1}(1\gamma _{n+2}\gamma _{n+1})_n+\gamma _{n+2}\gamma _{n+1}^3\gamma _n^2_{n1}$$
(26)
where $`_0=1`$, $`_1=12\gamma _0\gamma _1`$, $`_2=(1\gamma _1\gamma _2)^22\gamma _0\gamma _1(1\gamma _1\gamma _2)`$, and $`\gamma _n=\eta /[2(\delta n^2)]`$ for a solution with the period-$`2\pi `$ and $`_1=1\gamma _1^2`$, $`_2=(1\gamma _1\gamma _2)^2\gamma _1^2`$, $`_3=(1\gamma _1\gamma _2\gamma _2\gamma _3)^2\gamma _1^2(1\gamma _2\gamma _3)^2`$, and $`\gamma _n=2\eta /[4\delta (2n1)^2]`$ for a solution with the period-$`4\pi `$.

The stability diagram for Eq. (25) is shown in Fig. 5. The unstable regions, where at least one solution is unbounded, are shaded, and the stable regions, where all solutions are bounded, are not shaded. The boundary curves separating these regions are periodic solutions with period $`2\pi `$ (dashed lines) and $`4\pi `$ (solid lines). These curves are obtained by calculating $`_n=0`$ for $`n=200`$. The filled squares represent the values of $`(\overline{\delta }_m^T,\overline{\eta }_m^T)`$ satisfying the stability condition. Here we set $`\overline{\delta }_m^T=0`$ (i.e., $`\beta =0`$) since $`\overline{\delta }_m^T\beta ^2`$ and the terms of this order $`𝒪(\beta ^2)`$ have been neglected due to small $`\beta `$. For $`\overline{\delta }_m^T=0`$, the following values satisfy the stability condition: the values shown in Fig. 5 are $`0<\overline{\eta }_m^T<0.4540`$, $`\overline{\eta }_m^T`$3.7898, 10.6516, and 20.9637, and the values not shown in Fig. 5 are $`\overline{\eta }_m^T`$34.7142, 51.9022, 72.52784, 96.5910, 124.0918 $`\mathrm{}`$. For a large $`\overline{\eta }_m^T`$ satisfying the stability condition, we assume that $`\mathrm{\Omega }_m0`$ since $`\mathrm{\Omega }_m0`$ as $`\overline{\eta }_m^T\mathrm{}`$. In this case, the velocity for uniform motion is given by

$$\frac{\overline{\omega }_m}{k_x}\lambda _c\omega _p\left(\frac{1+2𝒮}{𝒜_m𝒞_m}\right)^{1/2},$$
(27)
indicating that the presence of the charging effect yields the mode dependent modification to the stability condition. For example, $`\overline{\omega }_1/k_x\lambda _c\omega _p`$ for $`m=1`$ (i.e., rectangular lattice) but $`\overline{\omega }_N/k_x\lambda _c\omega _p[(1+2𝒮)/(12𝒮)(1+4\alpha )]^{1/2}`$ for $`m=N`$ (i.e., triangular lattice). The velocity for the out-of-phase modes is reduced from the predicted value of the magnetic induction model (i.e., $`\alpha =0`$). Equation (27) indicates that moving Josephson vortices in a periodic array evolve from one stable mode to another as the vortex velocity increases. For example, as the vortex velocity exceeds $`\overline{\omega }_N/k_x`$, but less than $`\overline{\omega }_{N1}/k_x`$, the moving triangular lattice ($`m=N`$) becomes unstable and the $`m=N1`$ mode becomes stable.

## V Multiple sub-branching of supercurrent

In the resistive state, the supercurrent branch splits into multiple sub-branches as the bias current exceeds the Josephson current. Note that these supercurrent sub-branches differ from the observed multiple quasiparticle branches in the I-V data for LJJ stacks. This supercurrent sub-branching phenomenon, which appears clearly in the non-linear I-V regime, is attributed to the motion of Josephson vortices, but its origin is not understood clearly. Microwave induced voltage steps and geometric resonance are considered as other mechanisms, but we do not discuss them here. Instead, we argue that the splitting of the supercurrent branch is indeed due to the collective motion of Josephson vortices examining the low bias current regime where the I-V characteristics is linear. An analytic calculation is more tractable in this regime. Here, we illustrate qualitatively, rather than quantitatively, how the charging effect modifies the supercurrent sub-branches since the particle-hole imbalance effect neglected in this study may also need to be included for a quantitative comparison with the I-V data.

The current-voltage relations in the resistive state is obtained easily by noting that an AC voltage ripple with the Josephson frequency $`\omega _{\mathrm{},\mathrm{}1}`$, in addition to the DC voltage, appears across the junction when a bias current ($`J_B`$), greater than the critical current, is applied. This AC voltage ripple is due to the electron-pair tunneling current across the junction. Using the modified AC Josephson relation of Eq. (9), the time dependence of the phase difference between $`\mathrm{}`$ th and $`\mathrm{}1`$ th S layers can be written as

$$\phi _{\mathrm{},\mathrm{}1}(t)\phi _{\mathrm{},\mathrm{}1}(0)+\overline{\omega }_{\mathrm{},\mathrm{}1}t+\frac{\varphi _o}{2\pi }\frac{V_{\mathrm{},\mathrm{}1}^s}{\overline{\omega }_{\mathrm{},\mathrm{}1}}\mathrm{sin}\omega _{\mathrm{},\mathrm{}1}t$$
(28)
where $`\overline{\omega }_{\mathrm{},\mathrm{}1}=(2\pi /\varphi _o)[V_{\mathrm{},\mathrm{}1}+\alpha (V_{\mathrm{}+1,\mathrm{}}2V_{\mathrm{},\mathrm{}1}+V_{\mathrm{}1,\mathrm{}2})]`$ is the Josephson frequency, $`V_{\mathrm{},\mathrm{}1}`$ is the DC voltage (i.e., time averaged) across the superconductor layers $`\mathrm{}`$ and $`\mathrm{}1`$. $`V_{\mathrm{},\mathrm{}1}^s`$ is the amplitude of the AC voltage ripple. This time dependent phase difference of Eq. (28) yields a DC critical current response of

$$J_c\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}(t)=\mathrm{J}_1\left(\frac{\varphi _o}{2\pi }\frac{V_{\mathrm{},\mathrm{}1}^s}{\overline{\omega }_{\mathrm{},\mathrm{}1}}\right)\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}(0)$$
(29)
across the $`\mathrm{}`$ th and $`\mathrm{}1`$ th $`S`$ layers, indicating that the junction becomes resistive when $`J_B`$ exceeds the DC critical current. Here $`\mathrm{J}_1(x)`$ is the first order Bessel function of the first kind. Equation (29) indicates that the current $`J_{\mathrm{},\mathrm{}1}=J_BJ_c\mathrm{sin}\phi _{\mathrm{},\mathrm{}1}(t)`$ between two adjacent S layers is not uniform along $`z`$-direction, even though a uniform bias current is applied. Hence, in this resistive state, we may reduce Eq. (14) to

$$\left(\frac{1}{\omega _p^2}\frac{^2\mathrm{\Delta }_{\mathrm{}}}{t^2}+\frac{\beta }{\omega _p}\frac{\mathrm{\Delta }_{\mathrm{}}}{t}\right)+\frac{\alpha }{\omega _p^2}\frac{^2\mathrm{\Xi }_{\mathrm{}}}{t^2}=\frac{1}{J_c}\left[d^{}J_{\mathrm{},\mathrm{}1}+s(J_{\mathrm{}+1,\mathrm{}}+J_{\mathrm{}1,\mathrm{}2})\right].$$
(30)
Here, we neglected the spatial dependence (i.e., $`x`$ variation) of $`\phi _{\mathrm{},\mathrm{}1}`$ for simplicity. To explicitly express the I-V relation for each collective mode, we now rewrite Eq. (30) in Fourier space for the $`z`$-direction as

$$\frac{𝒜_m}{\omega _p^2}\frac{^2\phi _m}{t^2}+\frac{\beta }{\omega _p}\frac{\phi _m}{t}=\frac{J_m}{J_c}.$$
(31)
Note that the first and the second term on the left hand side of Eq. (31) represent the capacitive and the resistive contribution of the junction, respectively.

We now average Eq. (31) over time. Since AC Josephson tunneling leads to a small voltage oscillation about the DC voltage, a further simplification of Eq. (31) can be made. The modified AC Josephson relation of $`\phi _m/t=(2\pi /\varphi _o)V_m𝒜_m`$ indicates that the capacitive contribution vanishes when it is averaged over time (i.e., $`^2\phi _m/t^2V_m/t0`$). This simplification leads to the current-voltage relation of

$$V_m=\frac{\omega _p\varphi _o}{2\pi \beta }\frac{\overline{J}}{J_c}\frac{1}{𝒜_m}.$$
(32)
Here $`\mathrm{}`$ denotes the time average, and $`\overline{J}=J_mJ_BJ_c\mathrm{sin}\phi `$. Since the collective modes for $`m=1,2,\mathrm{}N`$ are identical to the modes for $`m=2N+1,2N,\mathrm{}N+2`$, respectively, the number of sub-branches is the same as the number of junctions (i.e., $`N`$) in the stack.

In Fig. 6, we plot the I-V relation of Eq. (32) for (a) $`\alpha =0.0`$ and (b) 0.1 to illustrate the effect of weak, but non-zero, charging effect. For clarity, we plot the curves for only the three collective modes (i.e., $`m=1`$, $`N/2`$, and $`N`$) corresponding to the modes shown in Fig. 2. These I-V curves reveal two interesting points. First, the supercurrent splits into $`N`$ sub-branches, each corresponding to the collective mode, when the LJJ stack are in the resistive state. The $`m=1`$ mode represents the high velocity mode (i.e., rectangular lattice), while the $`m=N`$ mode represents the low velocity mode (i.e., triangular lattice). Second, these $`N`$ sub-branches appear as a single curve when the charging effect is absent (i.e., $`\alpha =0`$, see Fig. 6(a)), but they spread out when this effect is present (i.e., $`\alpha 0`$, see Fig. 6(b)), suggesting that this can be used as another way to measure the charging effect. Since the width of this spread is related to the strength of the charging effect ($`\alpha `$), identification of each sub-branches is feasible at a low bias current. When $`\alpha `$ is small, as in BSCCO (i.e., $`\alpha 0.10.4`$), observing the branch splitting in the linear I-V regime may be difficult but is still possible. The main difficulty is in observing the high velocity branches (i.e., $`m𝒪(1)`$). To observe these branches, a magnetic field, stronger than $`BH_o`$, may be needed because of their stability conditions. Note that the appearance of these high velocity branches is expected when the interaction between the vortices is increased by the field, suggesting that a complete sub-branch structure may be more easily obtained from the I-V characteristics of a LJJ stack with increasing microwave irradiation power (i.e., AC magnetic fields). These results are consistent with the data exhibiting supercurrent branch splitting.

## VI summary and conclusion

In summary, we investigated the collective phase dynamics in the conventional LJJ stacks and in layered superconductors, using a theoretical model which accounts for both the magnetic induction effect and the charging effect. These two coupling mechanisms are equally important in the intrinsic LJJ (e.g. BSCCO) due to the atomic length thick S layers. We showed that the collective phase dynamics in an intrinsic LJJ stack is modified from those in a conventional LJJ stack in two important ways. (i) The dispersion of Josephson plasma mode for BSCCO is significantly changed from the Nb-Al/AlOx-Nb multilayers. Consequently, the Swihart velocity and the velocity of stable uniform motion for the out-of-phase collective modes in BSCCO increases and decreases, respectively, from the results of the magnetic induction model due to the presence of the charging effect. (ii) The supercurrent sub-branching in the resistive state is consistent with collective motion of Josephson vortices. The width of spread of these supercurrent sub-branches in the linear I-V regime depends on the strength of the charging effect. These results are consistent with the experimental data and illustrate the importance of accounting for the charging effect in BSCCO, even though its strength is weak ($`\alpha 0.10.4`$). They also suggest that our model is useful for understanding the experimental data for JPR, non-Josephson emission, and the I-V characteristics in the resistive state. Since many applications of intrinsic LJJ stacks as high frequency devices exploit collective dynamics of Josephson vortices, these results indicate that our model is useful for future technological applications involving intrinsic LJJ stacks.

ACKNOWLEDGMENTS

J.H.K. would like to thank H. J. Lee, W. Schwalm and D. H. Wu for helpful discussions. This work was supported in part by the UND Research Seed Money Program.

FIGURE CAPTIONS

Figure 1. A stack of LJJ is shown schematically as alternating superconducting (S) and insulating (I) layers with thickness $`d_S`$ and $`d_I`$, respectively. $`L_x`$ denotes the dimension in $`x`$direction. The magnetic field $`B`$ is applied in the plane of tunnel barriers and the bias current $`I_B`$ is applied along the vertical stack.

Figure 2. Mutual phase-locking of Josephson vortices (filled ovals) with the wave number $`k_m=m\pi /(N+1)a`$ is schematically illustrated. $`N`$ is the number of LJJ in the stack and $`a=d_S+d_I`$. For the in-phase mode ($`m=1`$), the Josephson vortices form a rectangular lattice, but for the out-of-phase mode ($`m=N`$), they form a triangular lattice. The dotted lines are a guide to eyes for phase-locking, and the arrows indicate the direction of propagation.

Figure 3. The dispersion relation for (a) the longitudinal plasma excitations at $`k_x=0`$ and (b) the transverse plasma excitations at $`k_m=0`$ are plotted to illustrate the difference between (a) our model and the magnetic induction model, and (b) our model and the charging effect model. Here $`\alpha =0.2`$ and $`𝒮=0.49999`$ are chosen.

Figure 4. The plasmon dispersion in the Meissner state (i.e., $`\phi _{\mathrm{},\mathrm{}1}^{(0)}=0`$) is plotted as functions of $`k_ma`$ and $`k_x\lambda _c`$ for the parameters corresponding to (a) the Nb-Al/AlOx-Nb multilayers ($`\alpha =0.0`$, $`𝒮=0.47`$) and (b) Bi2Sr2CaCu2O8+y ($`\alpha =0.2`$, $`𝒮=0.49999`$).

Figure 5. Stability diagram of Mathieu’s equation. The unstable regions and the stable regions are shaded and not shaded, respectively. The periodic solutions of period $`2\pi `$ (dashed lines) and $`4\pi `$ (solid lines) represent the boundary between the stable and unstable regions. The filled squares represent the values of ($`\overline{\delta }_m^T`$, $`\overline{\eta }_m^T`$) satisfying the stability condition.

Figure 6. The I-V curves for the supercurrent branch in the resistive state are plotted for (a) $`\alpha =0.0`$ and (b) 0.1 to illustrate the effect of nonzero $`\alpha `$ (i.e., charging effect). Here, only three curves corresponding to the collective modes shown in Fig. 2 are plotted for clarity.


format: text

text:
# Notes on Quasi-Homogeneous Functions in Thermodynamics

## I introduction

Quasi-homogeneous functions have been introduced in the framework of standard thermodynamics with the aim to studying scaling and universality near the critical point . A common synonymous of “quasi-homogeneous function” is “generalized homogeneous function” \[see e.g. \]. We wish to point out here that quasi-homogeneity can be an useful tool in the framework of standard thermodynamics, when one considers intensive variables as independent variables for the equilibrium thermodynamics description of a system. In fact, homogeneity for the fundamental equation in the entropy representation \[and in the energy representation\] is well-defined in terms of the standard Euler theorem for homogeneous functions . One simply defines the standard Euler operator (sometimes called also Liouville operator) and requires the entropy \[energy\] to be an homogeneous function of degree one. When the other thermodynamic potentials which are obtained from the entropy \[energy\] are taken into account by means of suitable Legendre transformations, then part of the independent variables are intensive . The thermodynamic potentials are still homogeneous of degree one in the extensive independent variables, but a different rescaling is appropriate for the independent variables. For example, let us consider the Gibbs potential $`G(T,p,N)`$ for a system which is described by means of three independent variables $`T,p,N`$. $`G`$ is homogeneous of degree one when the system is rescaled by $`\lambda `$, such a rescaling corresponding only to a rescaling $`N\lambda N`$, because $`T`$ and $`p`$ are intensive and remain unchanged under rescaling of the system. This is evident because, as it is well known, one has $`G=\mu (p,T)N`$, where $`\mu `$ is the chemical potential. Actually, one could also define $`G`$ as a quasi-homogeneous function of degree one with weights $`(0,0,1)`$. Then the behavior under scaling is better defined. A mathematical treatment of the same problem is found in Ref. . The approach we present here is characterized by the more general setting allowed by the technology of quasi-homogeneous functions; the sections on the Gibbs-Duhem equation and the on Pfaffian forms contain a further analysis of some formal aspects of standard thermodynamics.

## II quasi-homogeneous functions and thermodynamics

Given a set of real coordinates $`x^1,\mathrm{},x^n`$ and a set of weights $`𝜶(\alpha _1,\mathrm{},\alpha _n)\text{I}\text{R}^n`$, a function $`F(x^1,\mathrm{},x^n)`$ is quasi-homogeneous of degree $`r`$ and type $`𝜶`$  if, under dilatations by a scale factor $`\lambda >0`$ one finds

$$F(\lambda ^{\alpha _1}x^1,\mathrm{},\lambda ^{\alpha _1}x^n)=\lambda ^rF(x^1,\mathrm{},x^n).$$
(1)
A differentiable quasi-homogeneous function satisfies a generalized Euler identity:

$$DF=rF,$$
(2)
where $`D`$ is the Euler vector field

$$D\alpha _1x^1\frac{}{x^1}+\mathrm{}+\alpha _nx^n\frac{}{x^n}.$$
(3)
Notice that (2) is a necessary and sufficient condition for a differentiable function to be quasi-homogeneous . It is also interesting to define quasi-homogeneous Pfaffain forms. A Pfaffian form

$$\omega =\underset{i=1}{\overset{n}{}}\omega _i(x)dx^i$$
(4)
is quasi-homogeneous of degree $`r\text{I}\text{R}`$ if, under the scaling

$$x^1,\mathrm{},x^n\lambda ^{\alpha _1}x^1,\mathrm{},\lambda ^{\alpha _n}x^n$$
(5)
one finds

$$\omega \lambda ^r\omega .$$
(6)
This happens if and only if the degree of quasi-homogeneity deg$`(\omega _i(x))`$ of $`\omega _i(x)`$ is such that deg$`(\omega _i(x))=r\alpha _ii=1,\mathrm{},n`$. For a discussion about quasi-homogeneity and for further references, see .

### A quasi-homogeneous potentials in standard thermodynamics

Let us consider a thermodynamic potential $`R(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)`$ depending on $`k`$ intensive variables $`y^1,\mathrm{},y^k`$ and $`nk`$ extensive variables $`x^{k+1},\mathrm{},x^n`$. $`R`$ is required to be quasi-homogeneous of degree 1 and its type is

$$𝜶=(\underset{k}{\underset{}{0,\mathrm{},0}},\underset{nk}{\underset{}{1,\mathrm{},1}}).$$
(7)
Then, one has

$$R=\underset{i=k+1}{\overset{n}{}}x^i\frac{R}{x^i}.$$
(8)
This expression of the thermodynamic potentials is well-known, it is sometimes referred to as the identity satisfied by the potentials at fixed intensive variables . A treatment on a mathematical ground of the same topic is found in Ref. . It is evident that, in order to ensure that $`R`$ is a degree one quasi-homogeneous function, the intensive variables can be at most $`n1`$, in which case (cf. also the following section)

$$R=x^n\frac{R}{x^n}x^nr(y^1,\mathrm{},y^{n1}),$$
(9)
where $`r(y^1,\mathrm{},y^{n1})`$ is of degree zero.
We recall that, given the fundamental equation of thermodynamics in the energy \[entropy\] representation, one can obtain other fundamental equations by means of the Legendre transform . It is easy to show that:
the Legendre transform with respect to a variable of weight $`\alpha `$ of a quasi-homogeneous function of degree $`r`$ is a quasi-homogeneous function of degree $`r`$ with the weight $`\alpha `$ changed into the weight $`r\alpha `$ of the Legendre-conjugate variable (theorem 2 of ).
Moreover,
the partial derivative with respect to a variable of weight $`\alpha `$ of a quasi-homogeneous function $`R`$ of degree $`r`$ is a quasi-homogeneous function of degree $`r\alpha `$ having the same type as $`R`$ (theorem 1 of ). See also .
These results allow to justify easily the following examples.
For the free energy $`F(T,V,N)`$, one has $`F=UTS`$, thus $`F`$ is a quasi-homogeneous function of degree 1 and of weights $`(0,1,1)`$, and

$$F(T,V,N)=V\frac{F}{V}+N\frac{F}{N}.$$
(10)
Analogously,

$$S(T,V,N)=V\frac{S}{V}+N\frac{S}{N}.$$
(11)
\[In fact, $`S=F/T`$ and theorem 1 of  can be applied\]. Moreover, given $`S(T,p,N)`$, one has

$$S(T,p,N)=N\frac{S}{N}.$$
(12)
In concluding this section, we point out that the distinction between degree and weights of thermodynamic variables is somehow artificial, a degree becoming a weight if the thermodynamic variable is changed into an independent variable (e.g., the degree zero of the pressure becomes a weight zero when $`p`$ is an independent variable).

## III Gibbs-Duhem equations

Herein we take into account the Gibbs-Duhem equations. Cf. also . Let us define

$$R_i\frac{R}{x^i};R_a\frac{R}{y^a}$$
(13)
one has

$$dR=\underset{a=1}{\overset{k}{}}R_ady^a+\underset{i=k+1}{\overset{n}{}}R_idx^i.$$
(14)
On the other hand, one obtains from (8)

$$dR=\underset{i=k+1}{\overset{n}{}}R_idx^i+\underset{i=k+1}{\overset{n}{}}x^idR_i.$$
(15)
The GD equation is then

$$\underset{a=1}{\overset{k}{}}R_ady^a\underset{i=k+1}{\overset{n}{}}x^idR_i=0.$$
(16)
This equation is related with the quasy-homogeneity symmetry of the potential. Let us define the Euler operator

$$X\underset{i=k+1}{\overset{n}{}}x^i\frac{}{x^i}.$$
(17)
Let us also define a 1-form

$$\omega _R\underset{a=1}{\overset{k}{}}R_ady^a+\underset{i=k+1}{\overset{n}{}}R_idx^i$$
(18)
where $`R_a`$ are quasi-homogeneous functions of degree one $`XR_a=R_a`$ and the $`R_i`$ are quasi-homogeneous functions of degree zero $`XR_i=0`$. Then $`\omega _R`$ is a quasi-homogeneous 1-form of degree one, in the sense that it satisfies $`L_X\omega _R=\omega _R`$, where $`L_X`$ is the Lie derivative associated with $`X`$. One can also define a function

$$Ri_X\omega _R,$$
(19)
where $`i_X`$ is the standard contraction operator. As a consequence, one finds

$$dR=d(i_X\omega _R)=i_Xd\omega _R+L_X\omega _R=i_Xd\omega _R+\omega _R$$
(20)
If $`\omega _R`$ is a closed 1-form (and then, exact in the convex thermodynamic domain), then $`d\omega _R=0`$ and $`dR=\omega _R`$, i.e. $`R`$ is the potential associated with $`\omega _R`$. Notice also that in the latter case one finds

$$i_Xd\omega _R=0$$
(21)
which corresponds to the Gibbs-Duhem equation. In fact, one has

$$d\omega _R=\underset{a=1}{\overset{k}{}}dR_ady^a+\underset{i=k+1}{\overset{n}{}}dR_idx^i$$
(22)
and

$`i_Xd\omega _R`$ $`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}(i_XdR_a)dy^a{\displaystyle \underset{a=1}{\overset{k}{}}}dR_a(i_Xdy^a)`$ (23)
$`+`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}(i_XdR_i)dx^i{\displaystyle \underset{i=k+1}{\overset{n}{}}}dR_i(i_Xdx^i)`$ (24)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}R_ady^a{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^idR_i=0,`$ (25)
where $`i_XdR_a=XR_a=R_a`$, and $`i_XdR_i=XR_i=0`$.
The converse is also true, i.e., if (21) is satisfied then from (20) follows that $`\omega _R`$ is closed.
The GD equation is then satisfied because of the equality of the mixed second derivatives of $`R`$ (Schwartz theorem) and because of the quasi-homogeneous symmetry. In fact, by defining $`Q_{\alpha \beta }`$ the matrix of the second partial derivatives of $`R`$, one finds

$$\underset{i=k+1}{\overset{n}{}}x^idR_i=\underset{a=1}{\overset{k}{}}\underset{i=k+1}{\overset{n}{}}x^iQ_{ia}dy^a+\underset{j=k+1}{\overset{n}{}}\underset{i=k+1}{\overset{n}{}}x^iQ_{ij}dx^j.$$
(26)
Then the Gibbs-Duhem equation (16) is equivalent to

$`{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^iQ_{ia}dy^a`$ $`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}R_ady^a`$ (27)
$`{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^iQ_{ij}dx^j`$ $`=`$ $`0.`$ (28)
The former formula (27) is implemented if both Schwartz theorem and the quasi-homogeneous symmetry are implemented. In fact,

$`{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^iQ_{ia}dy^a`$ $`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i\left({\displaystyle \frac{}{x^i}}{\displaystyle \frac{}{y^a}}R\right)dy^a`$ (29)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \frac{}{y^a}}\left({\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \frac{}{x^i}}R\right)dy^a`$ (30)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \frac{R}{y^a}}dy^a={\displaystyle \underset{a=1}{\overset{k}{}}}R_ady^a.`$ (31)
Also (28) is implemented, in fact

$`{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^iQ_{ij}dx^j`$ $`=`$ $`{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i\left({\displaystyle \frac{}{x^i}}{\displaystyle \frac{}{x^j}}R\right)dx^j`$ (32)
$`=`$ $`{\displaystyle \underset{j=k+1}{\overset{n}{}}}\left({\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \frac{}{x^i}}{\displaystyle \frac{R}{x^j}}\right)dx^j,`$ (33)
and the latter is zero because $`R/x^i`$ are functions of degree zero for all $`i=k+1,\mathrm{},n`$.

## IV Pfaffian form $`\delta Q_{rev}`$

Let us consider the Pfaffian form $`\delta Q_{rev}`$ for a system described by $`(T,V,N)`$, where $`T`$ is the absolute temperature; one has

$$\delta Q_{rev}=C_{VN}(T)dT+a(T,V,N)dV+b(T,V,N)dN.$$
(34)
$`\delta Q_{rev}`$ has to be integrable, i.e., it satisfies $`\delta Q_{rev}d(\delta Q_{rev})=0`$, and it is known that $`T`$ is an integrating factor for $`\delta Q_{rev}`$, with

$$\frac{\delta Q_{rev}}{T}=dS.$$
(35)
Then, one finds that

$$\frac{\delta Q_{rev}}{T}=\frac{C_{VN}(T)}{T}dT+\frac{a(T,V,N)}{T}dV+\frac{b(T,V,N)}{T}dN$$
(36)
is exact and a potential is given by

$$S=\frac{a(T,V,N)}{T}V+\frac{b(T,V,N)}{T}N.$$
(37)
Notice that the quasi-homogeneity of degree one of $`S`$ is the tool allowing to obtain this result. It is “trivial” that $`S`$ is the potential associated with $`\delta Q_{rev}/T`$, it is less trivial that its “homogeneity” leads to (37). For a proof, see the appendix. $`\delta Q_{rev}`$ is quasi-homogeneous of degree one and weights $`(0,1,1)`$. From the theory of quasi-homogeneous integrable Pfaffian forms , it is known that an integrating factor is also given by

$$f=a(T,V,N)V+b(T,V,N)N.$$
(38)
The proof is found in Ref. . It is evident that

$$f=TS.$$
(39)
Analogously, one can consider $`(T,p,N)`$ as independent variables

$$\delta Q_{rev}=C_{pN}(T)dT+\eta (T,p,N)dp+\zeta (T,p,N)dN,$$
(40)
in which case

$$f=\zeta (T,p)N=TS.$$
(41)
## A potentials of exact quasi-homogeneous Pfaffian forms

We show that, if

$$\omega =\underset{a=1}{\overset{k}{}}B_a(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dy^a+\underset{i=k+1}{\overset{n}{}}B_i(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^i$$
(A1)
is a $`C^2`$ exact quasi-homogeneous Pfaffian form of degree one, with $`B_a,x^i`$ quasi-homogeneous of degree one and $`B_i,y^a`$ quasi-homogeneous of degree zero with respect to the Euler operator

$$Y=\underset{i=k+1}{\overset{n}{}}x^i\frac{}{x^i},$$
(A2)
then

$$P(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)\underset{i=k+1}{\overset{n}{}}B_i(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)x^i$$
(A3)
is a potential associated with $`\omega `$. In fact, let us consider

$`dP`$ $`=`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}B_i(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^i+{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^idB_i(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)`$ (A4)
$`=`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}B_i(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^i+{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \frac{B_i}{x^j}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^j`$ (A5)
$`+`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \frac{B_i}{y^a}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dy^a.`$ (A6)
The exactness of the Pfaffian form $`\omega `$ implies that $`d\omega =0`$ and, in particular

$`{\displaystyle \frac{B_i}{y^a}}`$ $`=`$ $`{\displaystyle \frac{B_a}{x^i}}a=1,\mathrm{},k;i=k+1,\mathrm{},n,`$ (A7)
$`{\displaystyle \frac{B_i}{x^j}}`$ $`=`$ $`{\displaystyle \frac{B_j}{x^i}}i,j=k+1,\mathrm{},n`$ (A8)
Then, one obtains

$`{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \frac{B_i}{y^a}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dy^a`$ $`=`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{a=1}{\overset{k}{}}}{\displaystyle \frac{B_a}{x^i}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dy^a`$ (A9)
$`=`$ $`{\displaystyle \underset{a=1}{\overset{k}{}}}\left({\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \frac{}{x^i}}B_a(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)\right)dy^a`$ (A10)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{k}{}}}B_a(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dy^a,`$ (A11)
because each $`B_a`$ is quasi-homogeneous of degree one. On the other hand, one has

$`{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \frac{B_i}{x^j}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^j`$ $`=`$ $`{\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \underset{j=k+1}{\overset{n}{}}}{\displaystyle \frac{B_j}{x^i}}(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)dx^j`$ (A12)
$`=`$ $`{\displaystyle \underset{j=k+1}{\overset{n}{}}}\left({\displaystyle \underset{i=k+1}{\overset{n}{}}}x^i{\displaystyle \frac{}{x^i}}B_j(y^1,\mathrm{},y^k,x^{k+1},\mathrm{},x^n)\right)dx^j`$ (A13)
$`=`$ $`0,`$ (A14)
because each $`B_i`$ is quasi-homogeneous of degree zero.


text:
# 1 The Feynman graphs that contribute to the vertex correction to the form factors 𝐹⁰⁰ and 𝐹²¹ and the wave-function correction to 𝐹⁰⁰. 𝑓 (𝑓') denotes the up (down) type fermion.

HEPHY-PUB 765/02

hep-ph/0210330

Radiative corrections to single Higgs boson
production in $`e^\mathbf{+}e^{\mathbf{}}`$ annihilation

H. Eberl 111Speaker, W. Majerotto, and V. C. Spanos

Contribution to SUSY02, 10th International Conference on Supersymmetry and Unification of Fundamental Interactions, 17–23 June 2002, DESY Hamburg, Germany.

## Abstract

For energies relevant to future linear colliders, $`\sqrt{s}500\mathrm{GeV}`$, the $`WW`$ fusion channel dominates the Higgs boson production cross section $`e^+e^{}\overline{\nu }\nu h^0`$. We have calculated the one-loop corrections to this process due to fermion and sfermion loops in the context of the MSSM. As a special case, the contribution of the fermion loops in the SM has also been studied. In general, the correction is negative and sizeable of the order of $`10\%`$, the bulk of it being due to fermion loops.

No Higgs boson could be detected so far. The four LEP experiments delivered lower bounds  for the Standard Model (SM) Higgs mass, $`m_h114\mathrm{GeV}`$ and the light CP even Higgs boson mass $`m_{h^0}88.3\mathrm{GeV}`$ of the Minimal Supersymmetric Standard Model (MSSM). In $`e^+e^{}`$ collisions, for energies $`200\mathrm{GeV}`$, the production of a single Higgs boson plus missing energy starts to be dominated by $`WW`$ fusion , that is $`e^+e^{}\overline{\nu }_e\nu _eWW\overline{\nu }_e\nu _eh^0`$, whereas the Higgsstrahlung process  $`e^+e^{}Zh^0\overline{\nu }\nu h^0`$ becomes less important. The rates for the $`ZZ`$ fusion are generally one order of magnitude smaller than those of the $`WW`$ channel.

At LHC, in $`pp`$ collisions, the gluon-gluon fusion mechanism provides the dominant contribution to Higgs production. Recently, it has been argued that also the channels $`WWh^0/H^0\tau \overline{\tau }`$ and $`WW`$, can serve as suitable search channels at LHC even for a Higgs boson mass of $`m_h120\mathrm{GeV}`$ . It was also shown  that the $`ep`$-option at LHC would offer the best opportunity to search for a Higgs boson in the mass range $`m_h<140\mathrm{GeV}`$ with $`WW`$ (and $`ZZ`$) fusion as the most important Higgs boson production mechanism there. At Tevatron, with $`p\overline{p}`$ collisions at $`2\mathrm{TeV}`$, the $`WW`$ fusion process plays a less important rôle at least for $`m_h180\mathrm{GeV}`$ .

This contribution follows Ref.  where the leading one-loop corrections to the $`WWh^0`$ vertex in the MSSM were calculated. Because of their Yukawa couplings, the fermion/sfermion loops are taken into account. Then the total cross section in $`e^+e^{}`$ annihilation was worked out for $`e^+e^{}\overline{\nu }_e\nu _eWW\overline{\nu }_e\nu _eh^0`$. We also included the Higgsstrahlung process $`e^+e^{}Zh^0\overline{\nu }\nu h^0`$ and the interference between these two mechanisms. Because the Higgsstrahlung process is much smaller in this range, we have neglected its radiative corrections. We have also discussed the SM case.

As for energies $`\sqrt{s}>500\mathrm{GeV}`$ the dominant channel $`e^+e^{}\overline{\nu }\nu h^0`$ is by far the $`WW`$ fusion, taking into account only (s)fermion loops, the renormalization of the five-point function simplifies to the renormalization of the $`WWh^0`$ vertex with off-shell $`W`$ bosons. The renormalization of the other two vertices in the process (e. g. the $`e^{}\nu _eW^+`$ coupling) is absorbed. For the renormalization procedure the on-shell scheme has been adopted.

The one-loop part of the $`WWh^0`$ coupling can be expressed in terms of all possible form factors

$$\left(\mathrm{\Delta }g_{WW}^{h^0}\right)^{\mu \nu }=F^{00}g^{\mu \nu }+F^{11}k_1^\mu k_1^\nu +F^{22}k_2^\mu k_2^\nu +F^{12}k_1^\mu k_2^\nu +F^{21}k_2^\mu k_1^\nu +iF^ϵϵ^{\mu \nu \rho \delta }k_{1\rho }k_{2\delta },$$
(1)
where $`k_{1,2}`$ denote the four-momenta of the off-shell $`W`$-bosons. The full analytic expressions of the renormalized form factor $`F^{00}`$ and of $`F^{21}`$, which are based on the evaluation of the Feynman graphs given in Fig. 1 with additional counter terms, can be found in Ref. .

The Feynman graphs for the one-loop corrected amplitude are drawn in Fig. 2. The tree-level part of the squared amplitude was already calculated in Ref. , and the expresions for the one-loop part, $`2\mathrm{}\left[\mathrm{\Delta }_W\left(_W^{\mathrm{tree}}\right)^{}\right]`$ and $`2\mathrm{}\left[\mathrm{\Delta }_W\left(_Z^{\mathrm{tree}}\right)^{}\right]`$, can be found in Ref. .

For the calculation of the cross section at tree-level, it is possible to perform some of the phase space integrations analytically and the rest of them numerically . However, including the one-loop correction terms, it is impossible to perform any of these integrations analytically .

Now let us discuss to the numerical results. The tree-level $`WWh^0`$ coupling for values of $`\mathrm{tan}\beta >5`$ as preferred by the LEP Higgs boson searches, mimics the SM one. For the calculation of the fermion/sfermion one-loop corrections to the $`WWh^0`$ vertex, the contribution of the third family of fermions/sfermions has been taken into account. This contribution turns out to be the dominant one, in comparison with the first two families corrections, due to the large values of the Yukawa couplings $`h_t`$ and $`h_b`$. The effect of the running of the coupling constants $`g`$ and $`g^{}`$ has been taken into account.

For simplicity, for all plots we have used $`A_t=A_b=A_\tau =A`$, $`\{m_{\stackrel{~}{U}},m_{\stackrel{~}{D}},m_{\stackrel{~}{L}},m_{\stackrel{~}{E}}\}=`$ $`\{\frac{9}{10},\frac{11}{10},1,1\}m_{\stackrel{~}{Q}}`$ and $`M_1=\frac{5}{3}M_2\mathrm{tan}^2\theta _W`$. The choice of a common trilinear coupling and the correlation between the soft sfermion masses are inspired by unification.

In Fig. 3 the parameters $`\mathrm{tan}\beta =10`$, $`\mu =100\mathrm{GeV}`$, $`A=500\mathrm{GeV}`$, $`m_{\stackrel{~}{Q}}=300\mathrm{GeV}`$, $`M_A=500\mathrm{GeV}`$, and $`M_2=400\mathrm{GeV}`$ are taken. Note that choosing different sets of parameters, the basic characteristics of these plots remain indifferent. In the left figure the various cross sections as a function of $`\sqrt{s}`$ for values up to $`1\mathrm{TeV}`$ are plotted. The dotted-dashed line represents the contribution from the $`WW`$ channel at tree-level alone, whereas the dotted line includes the Higgsstrahlung contribution as well. The dashed line comprises in addition the interference between the $`WW`$ channel and Higgsstrahlung. One can perceive that the size of this interference term is extremely small, and for this reason the difference between the dotted and dashed lines is rather minute. For $`\sqrt{s}500\mathrm{GeV}`$ the $`WW`$ fusion contribution dominates the total cross section for the Higgs production $`e^+e^{}\overline{\nu }\nu h^0`$. Actually, for $`\sqrt{s}800\mathrm{GeV}`$ the total tree-level cross section is due to $`WW`$ fusion. In the solid line we have taken into account the one-loop correction from the fermion/sfermion loops. The right figure of Fig. 3 shows the tree-level cross section $`\sigma _0`$ (dashed line) and the one-loop corrected cross section $`\sigma `$ (solid line) for energies up to $`3\mathrm{TeV}`$. Both figures show that the correction is always negative and in the order of $`10\%`$.

In Fig. 4 shows the relative correction $`\mathrm{\Delta }\sigma /\sigma _0`$ as a function of $`\sqrt{s}`$. The solid line corresponds to the set $`\mathrm{tan}\beta =10`$, $`\mu =100\mathrm{GeV}`$, $`A=500\mathrm{GeV}`$, $`m_{\stackrel{~}{Q}}=300\mathrm{GeV}`$, $`M_A=500\mathrm{GeV}`$, and $`M_2=400\mathrm{GeV}`$, whereas for the dashed line $`\mathrm{tan}\beta =40`$, $`\mu =300\mathrm{GeV}`$ and $`A=100\mathrm{GeV}`$ was taken, keeping the rest of parmeters unchanged. We see that the size of the one-loop correction to the Higgs production cross section becomes practically constant for $`\sqrt{s}>500\mathrm{GeV}`$ and weighs about $`15\%`$, almost independently of the choice of the SUSY parameters. This is due to the fact that the one-loop corrections are dominated by the fermion loops, and therefore the total correction is not very sensitive to the choice of the SUSY parameters.

In Fig. 5 we have plotted the cross section as a function of $`m_h`$ for the SM case, for $`\sqrt{s}=0.8\mathrm{TeV}`$ (red lines) and $`1\mathrm{TeV}`$ (black lines). The dashed lines correspond to the tree-level cross section for $`e^+e^{}\overline{\nu }\nu h^0`$, whereas the solid lines contains the one-loop correction stemming from the fermion loops. In addition, the couplings have been adjusted to the SM corresponding couplings. The plot exhibits the expected dependence of the cross section on $`m_h`$. What must be noticed is that especially for small Higgs boson masses $`200\mathrm{GeV}`$, the size of the fermion loops correction becomes important for the correct determination of the Higgs boson mass.

Finally, Fig. 6 exhibits the percentage of the sfermion loops to the total one-loop correction as a function of $`\mathrm{tan}\beta `$ (left) and $`\mu `$ (right), for two different values of $`\mu `$ and $`\mathrm{tan}\beta `$, respectively, as shown in the figure. In the left (right) figure $`A=100\mathrm{GeV}`$ ($`A=400\mathrm{GeV}`$). Further we have $`m_{\stackrel{~}{Q}}=300\mathrm{GeV}`$, $`M_A=500\mathrm{GeV}`$, and $`M_2=400\mathrm{GeV}`$, and $`\sqrt{s}`$ has been fixed to $`1\mathrm{TeV}`$. The grey area in the right figure is excluded due to the chargino mass bound. One can see that the maximum value of order 10% can be achieved for large values of $`\mu `$ and $`\mathrm{tan}\beta `$. There, due to the significant mixing in the stop and sbottom sector, the contribution of stops and sbottoms in the loops is enhanced. Even for such a parameter set the dominant correction, at least $`90\%`$ of the total correction, is due to the fermion loops.
An effective approximation of the one-loop corrections to the $`WWh^0`$ vertex can be found in , but it does not fully account for the whole effect.

In conclusion, we have calculated the fermion/sfermion loops corrections to the single Higgs boson production $`e^+e^{}\overline{\nu }\nu h^0`$ in the context of the MSSM and SM. They are supposed to be the dominant radiative corrections. For energies relevant to the future linear colliders, $`\sqrt{s}500\mathrm{GeV}`$, the $`WW`$ fusion channel dominates the cross section. In general, the correction due to fermion/sfermion loops is negative and yields a correction to the cross section of the order of $`10\%`$. The bulk of this correction stems from the fermion loops, and usually turns to be more than $`90\%`$ of the total correction. For the case of maximal mixing in the sfermion mass matrices, the contribution of the sfermion loops is enhanced, but nevertheless weighs less than $`10\%`$ of the total one-loop correction. As the correction is dominated by fermion loops and is rather independent of $`\sqrt{s}`$ for $`\sqrt{s}>500\mathrm{GeV}`$, we think that it can be approximated by a factor correction to the tree-level cross section. Such an approximation would be most useful for including initial state radiation (ISR) and beamstrahlung in an efficient way.

Acknowledgements
V. C. S. acknowledges support by a Marie Curie Fellowship of the EU programme IHP under contract HPMFCT-2000-00675. The authors acknowledge support from EU under the HPRN-CT-2000-00149 network programme and the “Fonds zur Förderung der wissenschaftlichen Forschung” of Austria, project No. P13139-PHY.


text:
# Zone edge focused two-phonon processes in He atom scattering from a simple prototype system: Xe(111).

## Abstract

We present a theoretical analysis of the multiphonon features appearing in the angular and energy resolved scattering spectra of low energy He atoms inelastically scattered from surface phonons. By applying the recently developed fully quantum multiphonon scattering formalism to the benchmark system Xe(111) we are able to reproduce with great accuracy the experimental data and also demonstrate how in the transition from a single- to a few- to a multiphonon scattering regime the one- and the many-phonon components of the scattering spectra evolve as the mean number of exchanged phonons increases with the increase of He atom incident energy. We show how the two-phonon processes may give rise to distinct peaks in the scattering spectrum as a result of the zone edge focusing effects and how the many-phonon features may still produce a structured background before the true multiphonon limit characterised by a Gaussian spectral shape is reached. These findings should prove useful in the interpretation of multiphonon He atom scattering spectra from other systems with similar surface vibrational properties.

Keywords: Atom-solid interactions, scattering; Energy dissipation; Phonons; Low index single crystal surfaces.

The energy and momentum transfer in gas-surface collisions proceeds under standard conditions dominantly in the multiphonon scattering regime. A convenient laboratory tool for studying these processes is inelastic scattering of thermal energy He beams from surfaces either in the ultra-high vacuum (UHV) environment, which enables to achieve and control the cleanliness of the studied surface on the atomic scale by means of the various surface science techniques, or in the wind tunnel experiments in which the energy transfer is directly obtainable from calorimetric measurements performed on samples exposed to the beam. However, the structure and cleanliness of the sample surface in the wind tunnel environment can be controlled and maintained only at a technical level.

The majority of UHV experiments utilising low energy He-atom-scattering time-of-flight (HAS-TOF) spectroscopy of surfaces have been performed with the aim to investigate and reveal the properties of surface phonons. To this end the experiments are carried out in the one-phonon scattering regime as only such experimental conditions allow clear fingerprinting of the surface phonon dispersion curves which serve as prerequisites for mode assignments and theoretical interpretations of the spectral peak intensities. The interpretation of HAS experiments in the multiphonon scattering regime is much more involved because of the need to deconvolute the various features in HAS-TOF spectra whose energies and intensities depend on phonon dispersions and scattering dynamics in a highly nontrivial fashion. Due to this complication very few experiments carried out in the multiphonon scattering regime have been complemented by adequate theoretical analyses, except in the case of high enough He beam energies for which the HAS-TOF spectra acquire a structureless Gaussian shape typical of the semiclassical scattering limit. A particularly complicated situation regarding the interpretation of measurements arises in the transition or cross-over from a single to a few-phonon scattering regime because of the difficulties encountered in theoretical assessments of the scattering events of this kind.

An exception to the general situation characterised by the difficulties to study transitions from a single- to a few- to a multiphonon scattering regime occurs for surfaces which exhibit nondispersive or Einstein-like phonons whose frequency $`\omega (𝐐)`$ as a function of a two-dimensional phonon wavevector $`𝐐`$ parallel to the surface remains constant over the surface Brillouin zone (SBZ), i.e. $`\omega (𝐐)=\omega _0`$. Typical examples of the systems exhibiting nondispersive surface phonons are monolayers of rare gas atoms adsorbed on low index crystal surfaces of metals and layers of molecular adsorbates like CO on metals. The HAS TOF spectra originating from single and multiphonon excitation of Einstein modes display a relatively simple structure consisting of a series of equidistant peaks located at multiples of Einstein phonon energy $`\mathrm{}\omega _0`$. The intensities of higher order peaks increase as the scattering regime changes from a single- to a multiphonon one but the actual transition between the two regimes depends on several parameters such as the incident energy and angle of the projectile atom, the magnitude of the mode frequency and the direction of its polarisation, surface temperature etc.

The quantity that best reflects the transitions between different scattering regimes is the exponent $`2W`$ of the Debye-Waller factor (DWF), $`\mathrm{exp}(2W)`$, that is calculated for a particular scattering spectrum. The value of $`2W=\overline{n}`$ is a function of the afore mentioned parameters and measures the mean number of phonons $`\overline{n}`$ exchanged in the course of collision. Hence, the transition from a single to a multiphonon scattering regime occurs for those values of the system and collision parameters for which $`2W1`$. In the case of nondispersive Einstein phonons this value of the Debye-Waller exponent (DWE) gives rise to maximum intensity of the first Einstein phonon peak in the scattering spectrum, whereas for $`2W2`$ the intensity of higher order peaks takes over. The recently developed quantum theory of multiphonon HAS from nondispersive phonons, which enables a detailed analysis of the properties of the DWF, has been successfully applied to interpret the TOF spectra of He atoms scattered from monolayers of Xe atoms on Cu substrates in which the dominant inelastic structure arises from excitation of Einstein modes. Many conclusions reached in these analyses can be conveniently extended also to the analyses of more complicated situations of inelastic HAS from dispersive phonons.

The problem of observation and pinpointing of few- and many-phonon processes in inelastic HAS from dispersive phonons is much more critical than in the case of Einstein modes and until recently no clearcut experimental evidence as to their resolution in the HAS TOF spectra was available. However, recent high resolution HAS studies of phonons characteristic of the (111) surface of fcc van der Waals crystal of xenon atoms condensed on Pt(111) substrate have revealed some distinct features in the TOF spectra which could not be attributed to the single phonon peaks over a structureless background. The dynamical matrix analysis of the Xe(111) surface enabled a clear assignment of the Rayleigh wave (RW), the longitudinal resonance (LR) and the shear horizontal (SH) phonons in the TOF spectra but not of the extra features which for incident He atom energies of $`E_i=10.5`$ meV and small parallel momentum transfer $`\mathrm{\Delta }𝐊`$ showed up as discernible maxima at energy loss $`\mathrm{\Delta }E4.2`$ meV.

To interpret the experimental results we have carried out fully quantum multiphonon calculations of the energy and parallel momentum resolved scattering spectra $`N_{𝐤_𝐢,T_s}(\mathrm{\Delta }E,\mathrm{\Delta }𝐊)`$ for experimentally fixed He atom initial momentum $`\mathrm{}𝐤_𝐢`$ and the substrate temperature $`T_s`$, using the exponentiated Born approximation (EBA) formalism. In this formalism the scattering spectrum as the function of the energy and parallel momentum transfer from the projectile to the phonons is given by a compact expression:

$$N_{𝐤_𝐢,T_s}^{EBA}(\mathrm{\Delta }E,\mathrm{\Delta }𝐊)=e^{2W}_{\mathrm{}}^{\mathrm{}}\frac{d\tau d^2𝓡}{(2\pi \mathrm{})^3}e^{\frac{i}{\mathrm{}}[(\mathrm{\Delta }E)\tau \mathrm{}(\mathrm{\Delta }𝐊)𝓡]}\mathrm{exp}[2W_{𝐤_𝐢,T_s}^{EBA}(𝓡,\tau )],$$
(1)
where the Debye-Waller exponent is $`2W=2W_{𝐤_𝐢,T_s}^{EBA}(0,0)`$ and the scattering function $`2W_{𝐤_𝐢,T_s}^{EBA}(𝓡,\tau )`$ is calculated in the EBA (for details see Eq. (12) of Ref. ). The values of $`\mathrm{\Delta }E`$ and $`\mathrm{\Delta }𝐊`$ appearing in expression (1) are not independent but are connected by the requirement of total energy and parallel momentum conservation in the collision, i.e. they are confined to the scan curve.

The calculated values of the Debye-Waller exponents pertaining to the measured spectra reported in Ref.  (c.f. Fig. 12 of Ref. ) indicated that these experiments were carried out in the two-phonon scattering regime. A typical result of the EBA calculations for the scattering spectrum of the He$``$Xe(111) collision system and comparison with the experimental data reported in Ref.  is shown in Fig. 1. Under these scattering conditions the observed features at $`\mathrm{\Delta }E=\pm `$ 4.2 meV, which can not be interpreted as arising from single phonon scattering processes, can be explained within the present model as a result of excitation of two RW phonons with wavevectors from nearly the opposite edges of the first SBZ so as that the total momentum transfer to phonons is small. A similar effect has also been observed in inelastic HAS from frustrated translation modes of c$`(2\times 2)`$ structure of CO on Cu(100). Following such a good agreement between the measured and theoretical results and taking into account the information on the scattering regime embodied in the corresponding DWE, we find a strong support for the above proposed assignments of the extra features in the He$``$Xe(111) TOF spectra. The goal of the present work is to further corroborate this interpretation through supplementing the previous calculations by more detailed theoretical analyses of the few-phonon scattering processes which may occur in HAS studies of this and other systems exhibiting similar surface vibrational properties.

In Fig. 2 we display four scattering spectra calculated for the collision system He$``$Xe(111) using the theoretical formalism outlined in Sec. V of Ref. . The parameter which is varied in calculations is the incident energy $`E_i`$ of the He atom whereas all other parameters are kept fixed (incident angle, surface temperature, etc). Since the maximum energies of the RW, LR and SH phonons that are localized at the Xe(111) surface are of the order of few meV only (see Fig. 11 of Ref. ), the probability of multiple excitation of surface phonons is non-negligible already for the lowest studied energy $`E_i=5.7`$ meV for which the scattering spectrum in the upper left-hand-side (LHS) panel of Fig. 2 has been calculated. This is corroborated by the magnitude of the corresponding Debye-Waller exponent $`2W=1.7`$ which implies that a part of the spectral weight has been shifted from the single to the multiphonon structure. The spectrum is dominated by the peaks signifying the excitation and annihilation of single RW and LR phonons in the scattering event but some structured multiphonon background is already discernible (dashed line) both under the one-phonon peaks and in the region between them. This multiphonon structure, and in the case studied this is a two-phonon structure, becomes more pronounced as the incoming energy is increased to $`E_i=10.4`$ meV for which $`2W=2.2`$, as shown in the upper right-hand-side (RHS) panel in Fig. 2. The peaks denoted by ”M” are of purely multiphonon origin but depending on their position they reflect kinematically different processes. The one at zero energy transfer $`\mathrm{\Delta }E=0`$ and denoted by M0 arises as a consequence of creation and annihilation of several phonons with zero total energy but nonzero total momentum balance. Note that the processes of this kind are possible only for finite substrate temperature and their effect in producing the high intensity of the no-loss or elastic line is particularly noticeable in the case of nondispersive phonons, as has been discussed in Refs.  and .

The origin of the peaks M2 and M-2 located at $`\mathrm{\Delta }E\pm 4.2`$ meV, respectively, is more subtle as it depends on a strong interplay between the characteristics of the phonon density of states typical of the studied system and the dynamics of the scattered particle. The RW phonons whose polarisation is dominantly perpendicular to the surface couple most strongly to the scattered He atoms, and this property holds in the single as well as in the multiphonon scattering regime. Since both the perpendicular polarization and the RW phonon density of states reach maximum at the SBZ boundary, it is expected that the largest contributions to the two-phonon processes will come from excitations of two Rayleigh waves with respective wavevectors $`𝐐_\mathrm{𝟏}`$ and $`𝐐_\mathrm{𝟐}`$ from the vicinity of the zone edges, i.e. $`𝐐_\mathrm{𝟏}=𝐐_\mathrm{𝟐}Q_{max}`$. On the other hand, as the one-phonon scattering matrix elements also depend on the projectile dynamics and produce maximum two-phonon intensity for minimum total momentum transfer, i.e. $`\mathrm{\Delta }𝐊=𝐐_\mathrm{𝟏}+𝐐_\mathrm{𝟐}=min`$, the combination of two requirements leads to the condition that maximum two-phonon scattering intensity occurs if the phonon wavevectors are pointed to the nearly opposite edges of the SBZ. Such constrained or ”zone edge focused” two phonon creation and annihilation processes give rise to the peaks in the scattering spectra that are located approximately at twice the value of the zone edge frequency of a single excited phonon. This is clearly visible in both upper panels of Fig. 2 in which the peaks M±2 indeed appear at $`\pm 4.2`$ meV, i.e. at twice the energy of RW phonon at the zone boundary (for dispersion curves of the RW phonon corresponding to the Xe(111) surface see Fig. 11 of Ref. ). The same effect occurs also in connection with excitation of LR and SH modes but is correspondingly smaller due to their much weaker coupling to the scattered He atom.

As the energy of the incident He atom is increased from $`E_i=10.4`$ meV to $`E_i=25`$ meV, the scattering regime changes from a two-phonon to a three-phonon one, as signified by the increase of the Debye-Waller exponent from the value $`2W=2.2`$ to the value $`2W=3.5`$. In the three-phonon inelastic He atom scattering processes the constraints on the wavevectors and energies of the emitted and absorbed phonons are much less severe and therefore the three-phonon structure appears much smoother than the two-phonon one. This is a general trend which holds for higher order phonon processes as well. The described situation is illustrated in the lower LHS panel of Fig. 2. Here the three-phonon structure, whose manifestation are the shoulders denoted by M±3 in the main body of the multiphonon spectrum, is the dominant one regarding the total spectral weight although the single- and two-phonon features still contribute significant components to the total intensity.

Lastly, as the incident projectile energy is increased to $`E_i=50`$ meV, the DWE reaches the value $`2W=7.2`$ and, accordingly, the calculated spectrum exhibits very little structure on top of a dominantly Gaussian-like shape typical of the multiphonon scattering regime. This is illustrated in the lower RHS panel in Fig. 2. The elastic peak, which is solely of the multiphonon origin, appears as a shoulder whereas the spectral maximum can be interpreted as a combination of many-phonon exchange processes resulting in the mean energy transfer of the order $`\mathrm{\Delta }\overline{E}=3`$ meV. Thus, the doubling of the DWE from 3.5 to 7.2 in going from the LHS to the RHS panel in Fig. 2 has resulted in a practically complete removal of the one- and few-phonon structures from the scattering spectra. This should not come as a surprise since the Gaussian spectral shape should always be expected from a distribution describing unitarized probabilities of uncorrelated scattering events in the limit of validity of the Stirling formula for approximating the factorials. Indeed, the last value of the DWE, viz. $`2W=7.2`$, lies at the lower boundary of this limit.

Finally, it may be of interest to investigate the validity of the Gaussian approximation for the multiphonon scattering formulae as given by Eqs. (73)-(75) of Ref.  and more accurately by Eqs. (291)-(293) of Ref. . To this end we have modelled the scattering spectrum for $`E_i=80`$ meV by the exact EBA expression from Eq. (10) of Ref.  and the Gaussian limit from Eq. (291) of Ref. . The results are shown in Fig. 3 and they illustrate that the Gaussian approximation is an excellent one under these scattering conditions and for the set of parameters characterising the He-Xe(111) interaction.

In summary, we have presented a detailed theoretical assessment of the energy and parallel momentum resolved scattering spectra of thermal energy He atoms inelastically scattered from surface phonons of a simple benchmark system Xe(111). The (111) surface of the van der Waals xenon crystal supports all three surface phonons, viz. the Rayleigh wave (RW), the longitudinal resonance (LR) and the shear horizontal (SH) mode, of which the RW due to its dominantly perpendicular polarization couples most strongly to the scattered He atoms. Hence, analogously as in the single phonon scattering regime, this mode should also give a dominant contribution to the multiple phonon exchange processes in the few- and multiphonon scattering regime. By varying the incoming energy of the scattered He atoms in expressions for the mean number of exchanged phonons $`\overline{n}`$ (equal to the Debye-Waller exponent $`2W`$) and for the inelastic scattering intensities we have numerically simulated the HAS-TOF spectra pertinent to the He$``$Xe(111) collision system in the various scattering regimes. We have shown by comparing experimental and theoretical results that in the transition from one phonon- to a two-phonon scattering regime the well resolved one-RW and one-LR phonon peaks represent the dominant structure both in the theoretical and measured spectra, but also that some extra peaks observed in experiments for incident He atom energy $`E_i=10.5`$ meV can be interpreted in terms of the ”zone edge focused” two-phonon processes. By increasing the incident He atom energy and thereby entering the regime of higher order phonon excitation processes ($`3\overline{n}6`$) we have shown that, contrary to the popular belief, the many-phonon contribution to the overall scattering spectrum may exhibit a relatively rich and complicated structure. Therefore, in this regime and before the multiphonon Gaussian limit is reached for $`\overline{n}7`$, the many-phonon contribution does not by any means represent a featureless background. Finally, the existence of a Gaussian limit is demonstrated by the calculations and comparisons of the exact EBA expression with the functional Gaussian limit of the same EBA multiphonon scattering spectrum for $`E_i=80`$ meV for which the mean number of exchanged phonons $`\overline{n}=10.9`$ is large enough to fullful the conditions for a true semiclassical multiphonon limit.

Although our findings and conclusions have been illustrated on a system that exhibits relatively simple vibrational properties, they should be of a more general validity. Hence they may serve as a guideline in pinpointing and detecting the peculiarities of the multiphonon structures also in other systems that have been studied by HAS-TOF spectroscopy in the various scattering regimes.

The authors would like to thank A.P. Graham and J.P. Toennies for communication and discussion of the experimental data on He atom scattering from Xe(111), and G. Witte and Ch. Wöll for useful comments on the properties of surface phonon density of states. This work has been supported in part by the Joint National Science Foundation grant JF 133 and the German-Croatian Collaboration Project-No.: KRO-007-97.