format: text

text:
# Hyperfine-mediated transitions between a Zeeman split doublet in GaAs quantum dots: The role of the internal field

## I Introduction

The electron spin states in bulk semiconductor and heterostructures have attracted much attention in recent years. Experiments indicate very long spin decoherence times and small transition rates between states of different spin kikkawa98:4313 ; ohno00:817 ; fujisawa01:R81304 . These promising results have motivated proposals for information processing based on electron spins in quantum dots, which might lead to a realization of a quantum computer loss98:120 ; burkard99:2070 .

A quantum dot is region where electrons are confined. The energy spectrum is discrete, due to the small size, and can display atomic-like properties ashoori96:413 ; tarucha96:3613 . Here we will consider quantum dots in GaAs-AlGaAs heterostructures. The main reasons for studying them are that relevant quantum dots are fabricated in such structures and Ga and As nuclei have a substantial hyperfine interactions with the conduction electrons. There are two main types of gate controlled dots in these systems, so-called vertical and lateral dotssohn97:105 . They differ in transverse confinement, which is approximately a triangular well and a square well for the lateral and vertical dots, respectively.

Manipulation the electron spin while maintaining phase coherence requires that it should be relatively well isolated from the environment. Coupling a quantum dot, or any closed quantum system, to its environment can cause decoherence and dissipation. One of the measures of the strength of the coupling to the environment are the transition rates, or inverse lifetimes, between the quantum dot states. In GaAs there are two main mechanisms that can cause finite lifetimes of spin states. These are the spin-orbit interaction and the hyperfine interaction with the surrounding nuclei. If a magnetic field is applied the change in Zeeman energy accompanying a spin-flip has to be absorbed by phonon emission. For two electron quantum dots, where the transitions are between triplet and singlet spin states, both spin-orbitkhaetskii00:12639 ; khaetskii00:470  and hyperfineerlingsson01:195306  mediated transitions have been studied. In both cases the transition rates are much smaller than the usual phonon rates, the spin-orbit rate being the higher except when the excited singlet and triplet states cross.erlingsson01:195306  For low magnetic fields, i.e. away from the single-triplet transition, the energy of the emitted phonon can be quite large and the transitions involve deformation phonons rather than piezoelectric ones.

If there is odd number of electrons in the dot the ground-state is usually a spin doublet so that energy change associated with the spin-flip is the electron Zeeman energy. Owing to the small $`g`$-factor in GaAs this energy is rather small compared to the level spacing and the dominating phonon mechanism is due to piezoelectric phonons. Recently spin-orbit mediated spin-flip transitions between Zeeman levels was investigatedkhaetskii01:27 . Due to Kramer’s degeneracy the transition amplitude for the spin-flip is proportional to the Zeeman splitting. This results in a spin-orbit spin-flip rate proportional to the fifth power of the Zeeman splitting.

In this publication we consider hyperfine mediated transitions between Zeeman split levels in a quantum dot. The transition amplitude remains finite at zero external magnetic field, giving a spin-flip rate (Eq. (27) ) that is proportional to the cube of the Zeeman splitting. The cause of this is an internal magnetic field due to the hyperfine interaction. We consider the important concept of internal magnetic field in some detail. Since the parameters of hyperfine interaction between conduction band electrons and underlying nuclei in GaAs have been extensively investigated paget77:5780 ; meier84:381 , including the Overhauser effect and spin-relaxation in GaAs/AlGaAs heterostructures dobers88:1650 ; berg90:2563 ; vagner82:635 ; barret95:5112 ; tycko95:1460  , we are able now to predict the typical time scale for this process in particular quantum dot geometries.

Upon completion of this work we learned about recent results of Khaetskii, Loss and Glazman khaetskii02:1 . They consider essentially the same situation and model and obtain electron spin decoherence without considering any mechanism of dissipation. This is in clear distinction from the present result for spin-flip rate that requires a mechanism of dissipation, i.e. phonons.

The rest of the paper is organized as follows: in section II the model used is introduced in addition to the basic assumptions and approximations used, section III deals with the internal magnetic field due to the hyperfine interaction and section IV contains the derivation of the transition rates. Finally, in section V the results are discussed.

## II Model and assumptions

We consider a quantum dot embedded in a AlGaAs/GaAs heterostructure. The details of the quantum dot eigenstates are not important for now, it suffices to say that the energy spectrum is discrete and the wave functions are localized in space. The spatial extension of the wavefunction in the lateral and transverse direction (growth direction) are denoted with $`\mathrm{}`$ and $`z_0`$, respectively. Quantum dots in these heterostructures are formed at a GaAs/AlGaAs interface, where the confining potential is very strong so that $`\mathrm{}z_0`$. We define the volume occupied by an electron as $`V_{\text{QD}}=\pi \mathrm{}^2z_0`$. The Hamiltonian of the quantum dot can be written on the form

$$H_0=\underset{l}{}\mathbf{\left(}\epsilon _l+g\mu _B𝑩\widehat{𝑺}\mathbf{\right)}|ll|,$$
(1)
where $`\epsilon _l`$ are the eigenenergies which depend on the structure of the confining potential and the applied magnetic field. The magnetic field also couples to the electron spin via the Zeeman term.

Since the $`\mathrm{\Gamma }`$ point of the conduction band in GaAs is mainly composed of $`s`$ orbitals the dipole interaction with the nuclei vanishes and the hyperfine interaction can be described by the usual contact term

$$H_{\text{HF}}=A\widehat{𝑺}\underset{k}{}\widehat{𝑰}_k\delta (𝒓𝑹_k),$$
(2)
where $`\widehat{𝑺}`$ ($`\widehat{𝑰}_k`$) and $`𝒓`$ ($`𝑹_k`$) denotes, respectively, the spin and position of the electron ($`k`$th nuclei). The delta function indicates the point-like nature of the contact interaction will result in a position dependant coupling. The coupling constant $`A`$ has the dimension Volume$`\times `$Energy. To get a notion of the related energy scale, it is straightforward to relate $`A`$ to the energy splitting of the doublet for a fully polarized nuclear system,

$$E_n=AC_nI,$$
(3)
$`C_n`$ being density of nuclei and $`I`$ is the spin of a nucleus. In GaAs this energy is $`E_n0.135`$ meV, which corresponds to a magnetic field of about 5 T.dobers88:1650  For a given quantum dot geometry the number of nuclei occupying the dot is defined as $`N_{\text{QD}}=C_nV_{\text{QD}}`$. Since the Ga and As nuclei have the same spin and their coupling constants are comparable, we will assume that all the nuclear sites are characterized by the same hyperfine coupling $`A`$, and the $`C_na_0^3`$ where $`a_0`$ is the lattice spacing. For realistic quantum dots $`N_{\text{QD}}10^410^6`$, and it is therefore an important big parameter in the problem.

The coupling between the electron and the phonon bath is represented with

$$H_{\text{ph}}=\underset{𝒒,\nu }{}\alpha _\nu (𝒒)(b_{𝒒,\nu }e^{i𝒒𝒓}+b_{𝒒,\nu }^{}e^{i𝒒𝒓}),$$
(4)
where $`b_{𝒒\nu }^{}`$ and $`b_{𝒒\nu }`$ are creation and annihilation operator for the phonon mode with wave-vector $`𝒒`$ on brach $`\nu `$. In GaAs there are two different coupling mechanisms, deformation ones and piezoelectric ones. For transitions between Zeeman split levels in GaAs, i.e. low energy emission, the most effective phonon mechanism is due to piezoelectric phonons. We will assume that the heterostructure is grown in the  direction. This is the case for almost all dots and it imposes important symmetry relations on the coupling coefficient. The square of the coupling coefficient for the piezoelectric phonons is then given by (see Ref. price81:217, )

$$\alpha _\nu ^2(𝒒)=\frac{(eh_{14})^2\mathrm{}}{2\rho c_\nu Vq}A_\nu (\theta )$$
(5)
where $`(eh_{14})`$ is the piezoelectric coefficient, $`\rho `$ is the mass density, $`c_\nu `$ is the speed of sound of branch $`\nu `$, $`V`$ is normalization volume and we have defined $`𝒒=q(\mathrm{cos}\varphi \mathrm{sin}\theta ,\mathrm{sin}\varphi \mathrm{sin}\theta ,\mathrm{cos}\theta )`$ and the $`A_\nu `$’s are the so-called anisotropy functions, see appendix B.

## III Internal magnetic field

In this section we will introduce the concept of the effective magnetic field, the internal field, acting on the electron due to the hyperfine interaction. This internal field is a semiclassical approximation to the nuclear system, this approximation being valid in the limit of large number of nuclei, $`N_{\text{QD}}1`$. If the nuclei are noticeably polarized, this field almost coincides with the Overhauser field that represents the average nuclear polarization. It is important that the internal field persists even at zero polarization giving rise to doublet splitting of the order $`E_nN_{\text{QD}}^{1/2}`$.

First we write the Hamiltonian in Eq. (2) in the basis of the electron orbital states and present it as a sum of two terms

$$H_{\text{HF}}=H_{\text{HF}}^0+V_{\text{HF}},$$
(6)
where the terms are defined as

$`H_{\text{HF}}^0`$ $`=`$ $`A{\displaystyle \underset{l}{}}|ll|\widehat{𝑺}{\displaystyle \underset{k}{}}|𝑹_k|l|^2\widehat{𝑰}_k`$ (7)
$`V_{\text{HF}}`$ $`=`$ $`{\displaystyle \underset{ll^{}}{}}|ll^{}|\widehat{𝑺}{\displaystyle \underset{k}{}}l|𝑹_k𝑹_k|l^{}\widehat{𝑰}_k.`$ (8)
By definition $`H_{\text{HF}}^0`$ does not couple different orbital levels. By combining Eqs. (1) and (7) one obtains the following Hamiltonian

$$H_0=\underset{l}{}\mathbf{\left(}\epsilon _l+g\mu _B\widehat{𝓑}_l\widehat{𝑺}\mathbf{\right)}|ll|.$$
(9)
We will regard the mixing term $`V_{\text{HF}}`$ as a perturbation to $`H_0`$. The justification for this is that the typical fluctuations of the electron energy due to the hyperfine interaction are much smaller than the orbital energy separation, $`\mathrm{}\mathrm{\Omega }`$.

We now concentrate on $`H_0`$ and formulate semiclassical description of it. For this we consider the operator of the orbitally dependent effective magnetic field;

$$\widehat{𝓑}_l=𝑩+\frac{1}{g\mu _B}\widehat{𝑲}_l$$
(10)
where

$$\widehat{𝑲}_l=\frac{E_n}{IC_n}\underset{k}{}|𝑹_k|l|^2\widehat{𝑰}_k.$$
(11)
Our goal is to replace the operator $`\widehat{𝑲}_l`$ with a classical field. To prove the replacement is reasonable we calculate the average of the square for a given unpolarized nuclear state $`|𝝁`$

$$𝑲_l^2=𝝁|\widehat{𝑲}_l^2|𝝁=\frac{E_n^2}{N_{\text{QD}}}\rho _l^{[2]}(I+1)/I$$
(12)
where $`\rho _l^{[n]}=𝑑𝒓V_{\text{QD}}^{n1}|𝒓|l|^{2n}`$ is a dimensionless constant depending on the wave function of the orbital $`l`$. We cannot simply replace $`\widehat{𝑲}_l`$ by its eigenvalues: as in the case of the usual spin algebra, different components of $`\widehat{𝑲}_l`$ do not commute. In addition its square does not commute with individual components, $`[\widehat{𝑲}_l^2,K_l^\alpha ]0`$. To estimate fluctuations of $`\widehat{𝑲}_l`$ we calculate the uncertainty relations between its components

$`\mathrm{\Delta }K_l^\alpha \mathrm{\Delta }K_l^\beta {\displaystyle \frac{E_n^2}{N_{\text{QD}}^{3/2}}}=𝑲_l^2{\displaystyle \frac{1}{N_{\text{QD}}^{1/2}}}.`$ (13)
Since $`N_{\text{QD}}^{1/2}1`$ we have proved that the quantum fluctuations in $`\widehat{𝑲}_l`$ are much smaller than its typical amplitude. The semiclassical picture introduced above is only valid for high temperatures, $`kTE_nN_{\text{QD}}^1`$, where there are many states available to the nuclear system and the typical length is proportional to $`N_{\text{QD}}^{1/2}`$. For temperatures below $`E_nN_{\text{QD}}^1`$ the nuclear system will predominantly be in the ground state and the classical picture breaks down. This is similar to the quantum mechanical description of a particle moving in a potential. At zero temperature it will be localized in some potential minimum and quantum mechanics will dominate. At sufficiently high temperatures the particle occupies higher energy states and its motion is well described by classical mechanics. Having established this we can replace the average over the density matrix of the nuclei by the average over a classical field $`𝑲_l`$, note that it has no ‘hat’, whose values are Gaussian distributed

$$P(𝑲_l)=\left(\frac{3}{2\pi 𝑲_l^2}\right)^{3/2}\mathrm{exp}\mathbf{\left(}\frac{3(𝑲_l𝑲_l^{(0)})^2}{2𝑲_l^2}\mathbf{\right)}.$$
(14)
Here $`𝑲_l^{(0)}`$ is the average, or Overhauser, field assuming small degree of polarization. The effective magnetic field acting on the electron is $`𝓑_𝒍=_l𝒏`$ where $`𝒏`$ is the unit vector along the total field for a given $`𝑲_l`$, see Fig. 1. For this configuration the spin eigenfunctions are $`|𝒏_\pm `$ corresponding to eigenvalues

$$𝒏\widehat{𝑺}|𝒏_\pm =\pm \frac{1}{2}|𝒏_\pm $$
(15)
The effective Zeeman Hamiltonian is then

$$H_\text{Z}=g\mu _B𝓑_l\widehat{𝑺}$$
(16)
in this given internal field configuration. The spectrum of $`H_0`$ thus consists of many doublets distinguished by the value of $`𝓑_l`$. The magnitude of the effective field $`_l`$ determines the Zeeman splitting of the each doublet

$`\mathrm{\Delta }_l`$ $`=`$ $`g\mu _B_l`$ (17)
$`=`$ $`(E_B^2+K_l^2+2E_BK_l\mathrm{cos}\theta )^{1/2}`$
where $`E_B=g\mu _BB`$ the external field Zeeman energy.

We conclude this section with two remarks concerning time and energy scales. First, since the dynamics of $`𝑲_l`$ is essentially the dynamics of the nuclear system it changes on a time scale of the nuclear relaxation, which can be quite long.berg90:2563  For electron dynamics on shorter timescales $`𝑲_l`$ plays the role of a ‘frozen disorder’. At longer time scales self-averaging over all configurations of $`𝑲_l`$ takes place. Second, the RMS value of $`𝑲_l`$ is approximately $`5`$$`\mu `$eV which corresponds to a 100 Gauss field (for a value $`N_{\text{QD}}10^5`$).

## IV Transition Rate

We concentrate on the transitions between the doublet components in Eq. (15). Assuming that the higher energy doublet state is initially occupied, we will calculate the transition rate to the lower one. The transition must be accompanied by energy dissipation equal to $`\mathrm{\Delta }_0`$. This energy cannot be absorbed be the nuclear system.kim94:16777 ; erlingsson01:195306  So we need an external mechanism of energy dissipation. The most effective one is known to be phonons. However the phonons alone cannot change the electron spin so we need a mechanism which mixes spin and orbital degrees of freedom, that is $`V_{\text{HF}}`$ from Eq. (8). Thus the transition amplitude is proportional to both $`V_{\text{HF}}`$ and the electron phonon coupling $`H_{\text{ph}}`$.

Here we assume that the electron is in the lowest orbital state $`|0`$ since the phonon mechanism will bring the electron to this state from any higher orbital on timescales much smaller than those related to transtions between the doublet components. Thus we consider an initial state of the entire system $`|\text{i}=|0,𝒏_{};𝝁;𝑵`$ which is a product state of the electron, nuclear $`|𝝁`$ and phonon $`|𝑵`$ systems and the final state $`|\text{f}=|0,𝒏_+;𝝁^{};𝑵^{}`$. The transition amplitude between them, in second order perturbation theory, reads

$`𝒯`$ $`=`$ $`{\displaystyle \underset{l0}{}}({\displaystyle \frac{0,𝒏_+;𝝁^{}|V_{\text{HF}}|l,𝒏_{};𝝁l;𝑵^{}|H_{\text{ph}}|0;𝑵}{(\epsilon _0\epsilon _l)+E_B}}`$ (18)
$`+{\displaystyle \frac{0;𝑵^{}|H_{\text{ph}}|l;𝑵l,𝒏_+;𝝁^{}|V_{\text{HF}}|0,𝒏_{};𝝁}{(\epsilon _0\epsilon _l)E_B}}).`$
The summation is over virtual states involving higher orbitals and the denominators in Eq. (18) contains the energy differences between different orbital states. The internal field depends on the orbital state, resulting in a rather complicated expression. Albeit the energy related to the internal field is much smaller than the orbital separation so we can safely replace the Zeeman splitting with $`E_B`$. The reason for that is that only at high external fields where $`\mathrm{\Delta }_0E_B`$ will the effects of the Zeeman splitting be appreciable in the denominator. The internal field also appears in the phonon rate since it determines the electron energy difference between the initial and final state. Since $`H_{\text{ph}}`$ does not connect different nuclear states and conversely $`V_{\text{HF}}`$ does not mix different phonon states the sums over intermediate phonon and nuclear states reduce to a single term. From this transition amplitude the transition rate is obtained via Fermis golden rule

$$\stackrel{~}{\mathrm{\Gamma }}_{\text{sf}}=\frac{2\pi }{\mathrm{}}\underset{𝑵^{}}{}\underset{𝝁^{}}{}|𝒯|^2\delta (E_\text{i}E_\text{f}),$$
(19)
where $`E_\text{i}E_\text{f}`$ is the energy difference between the initial and final state of the combined systems. Substituting Eq. (18) into Eq. (19) we get the following relation for the spin-flip rate

$`\stackrel{~}{\mathrm{\Gamma }}_{\text{sf}}`$ $`=`$ $`{\displaystyle \underset{l,l^{}0}{}}({\displaystyle \frac{l,𝒏_{};𝝁|V_{\text{HF}}|0,𝒏_+0,𝒏_+|V_{\text{HF}}|l^{},𝒏_{};𝝁}{((\epsilon _0\epsilon _l)+E_B)((\epsilon _0\epsilon _l^{})+E_B)}}{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝑵^{}}{}}0;𝑵|H_{\text{ph}}|l;𝑵^{}l^{};𝑵^{}|H_{\text{ph}}|0;𝑵\delta (E_\text{i}E_\text{f})`$ (20)
$`+{\displaystyle \frac{l,𝒏_{};𝝁|V_{\text{HF}}|0,𝒏_+l^{},𝒏_+|V_{\text{HF}}|0,𝒏_{};𝝁}{((\epsilon _0\epsilon _l)+E_B)((\epsilon _0\epsilon _l^{})E_B)}}{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝑵^{}}{}}0;𝑵|H_{\text{ph}}|l;𝑵^{}0;𝑵^{}|H_{\text{ph}}|l^{};𝑵\delta (E_\text{i}E_\text{f})`$
$`+{\displaystyle \frac{0,𝒏_{};𝝁|V_{\text{HF}}|l,𝒏_+0,𝒏_+|V_{\text{HF}}|l^{},𝒏_{};𝝁}{((\epsilon _0\epsilon _l)E_B)((\epsilon _0\epsilon _l^{})+E_B)}}{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝑵^{}}{}}l;𝑵|H_{\text{ph}}|0;𝑵^{}l^{};𝑵^{}|H_{\text{ph}}|0;𝑵\delta (E_\text{i}E_\text{f})`$
$`+{\displaystyle \frac{0,𝒏_{};𝝁|V_{\text{HF}}|l,𝒏_+l^{},𝒏_+|V_{\text{HF}}|0,𝒏_{};𝝁}{((\epsilon _0\epsilon _l)E_B)((\epsilon _0\epsilon _l^{})E_B)}}{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝑵^{}}{}}l;𝑵|H_{\text{ph}}|0;𝑵^{}0;𝑵^{}|H_{\text{ph}}|l^{};𝑵\delta (E_\text{i}E_\text{f})).`$
The spin-flip rate depends on the initial state of the nuclear system $`|𝝁`$. This poses the problem of how to deal with the nuclear state $`|𝝁`$, since we already demoted all the spin operators to a collective classical variable. A conceptually simple solution lies in the fact that when Eq. (20) and (8) are considered together one sees that the rate is a sum over all pairs of nuclei in the system. Focusing on a given pair of nuclei $`k`$ and $`k^{}`$, all the other nuclei are unchanged when the electron spin is ‘scattered’ on by this pair. By simply redefining the classical field such that it is composed of all nuclei except this given pair we can circumvent the problem. This procedure will not change our previous result regarding the properties of $`𝑲_l`$ and also defining $`|𝝁=|m_k|m_k^{}`$ makes it straightforward to work with the nuclear states in Eq. (20).

Although the transition rate can be very slow, the typical duration of a transition event is set by energy uncertainty, $`\mathrm{}/\mathrm{\Delta }_0`$. This is much shorter than the typical time for nuclear spin relaxation so that the nuclear system is frozen in a given configuration of $`𝑲_0`$. In this case taking an average over $`𝑲_0`$ is not well motivated since a given configuration of $`𝑲_0`$ is fixed. For now we will postpone the average over the classical field. After performing a thermal average over $`|𝝁`$ in the transition rate we obtain the following

$`\mathrm{\Gamma }_{\text{sf}}`$ $`=`$ $`G_{\text{corr}}({\displaystyle \underset{l0}{}}\{{\displaystyle \frac{2a_{ll}\gamma _{ll}}{\delta \epsilon _l^2}}(1+3{\displaystyle \frac{E_B^2}{\delta \epsilon _l^2}})+{\displaystyle \frac{2\mathrm{}\{\stackrel{~}{a}_{ll}\stackrel{~}{\gamma }_{ll}\}}{\delta \epsilon _l^2}}(1+{\displaystyle \frac{E_B^2}{\delta \epsilon _l^2}})\}`$ (21)
$`+{\displaystyle \underset{l<l^{}0}{}}\{{\displaystyle \frac{4\mathrm{}\{a_{ll^{}}\gamma _{ll^{}}\}}{\delta \epsilon _l\delta \epsilon _l^{}}}\mathbf{(}1+\mathbf{(}{\displaystyle \frac{(\delta \epsilon _l+\delta \epsilon _l^{})^2}{\delta \epsilon _l\delta \epsilon _l^{}}}1\mathbf{)}{\displaystyle \frac{E_B^2}{\delta \epsilon _l\delta \epsilon _l^{}}}\mathbf{)}+{\displaystyle \frac{4\mathrm{}\{\stackrel{~}{a}_{ll^{}}\stackrel{~}{\gamma }_{ll^{}}\}}{\delta \epsilon _l\delta \epsilon _l^{}}}\mathbf{(}1+\mathbf{(}{\displaystyle \frac{(\delta \epsilon _l\delta \epsilon _l^{})^2}{\delta \epsilon _l\delta \epsilon _l^{}}}+1\mathbf{)}{\displaystyle \frac{E_B^2}{\delta \epsilon _l\delta \epsilon _l^{}}}\mathbf{)}\})`$
where $`\delta \epsilon _l=\epsilon _0\epsilon _l`$. Here we have introduced the one-site nuclear spin correlation function $`G_{\text{corr}}`$, which is evaluated in appendix A. The denominator has been expanded to second order in the Zeeman splitting. The parameters $`a_{ll^{}},\stackrel{~}{a}_{ll^{}}`$ are related to the $`V_{\text{HF}}`$ matrix elements

$`a_{ll^{}}`$ $`=`$ $`A^2C_n{\displaystyle d^3𝑹_k\mathrm{\Psi }_l^{}(𝑹_k)|\mathrm{\Psi }_0(𝑹_k)|^2\mathrm{\Psi }_l^{}(𝑹_k)}`$ (22)
$`\stackrel{~}{a}_{ll^{}}`$ $`=`$ $`A^2C_n{\displaystyle d^3𝑹_k\mathrm{\Psi }_l^{}(𝑹_k)\mathrm{\Psi }_l^{}^{}(𝑹_k)\mathrm{\Psi }_0(𝑹_k)^2}`$ (23)
and $`\gamma _{ll^{}},\stackrel{~}{\gamma }_{ll^{}}`$ are generalized phonon transition rates

$`\gamma _{ll^{}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝒒\nu }{}}{\displaystyle \frac{\alpha _\nu ^2(𝒒)[e^{i𝒒𝒓}]_{0,l}[e^{i𝒒𝒓}]_{l^{},0}}{e^{\beta \mathrm{}\omega _{𝒒,\nu }}1}}\delta (\mathrm{}\omega _{𝒒\nu }\mathrm{\Delta }_0)`$ (24)
$`\stackrel{~}{\gamma }_{ll^{}}`$ $`=`$ $`{\displaystyle \frac{2\pi }{\mathrm{}}}{\displaystyle \underset{𝒒\nu }{}}{\displaystyle \frac{\alpha _\nu ^2(𝒒)[e^{i𝒒𝒓}]_{0,l}[e^{i𝒒𝒓}]_{0,l^{}}}{e^{\beta \mathrm{}\omega _{𝒒,\nu }}1}}\delta (\mathrm{}\omega _{𝒒\nu }\mathrm{\Delta }_0).`$ (25)
Here we have only included the emission process since we assume the the spin is initially in the higher energy state.

Until now we have considered a general quantum dot and the rate in Eq. (21) is valid for any quantum dot. To proceed further we will specify the confining potential to be parabolic in the lateral direction and in the transverse one a triangular well potential is chosen. The wave function is $`𝒓|l\chi _0(z)\psi _{n,M}(r,\theta )`$ where $`n,M`$ denote the orbital and angular momentum quantum numbers respectively of the Darwin-Fock solution and $`\chi _0(z)`$ is the wave function in the transverse direction. The square of the lateral confining length is $`\mathrm{}^2=\mathrm{}^2/m^{}\mathrm{}\mathrm{\Omega }`$, where $`m^{}`$ is the electron effective mass and $`\mathrm{\Omega }=(\mathrm{\Omega }_0^2+(\omega _c/2)^2)^{1/2}`$ is the effective confining frequency. The confining energy is $`\mathrm{}\mathrm{\Omega }_0`$ and $`\omega _c=eB/m^{}`$ is the cyclotron frequency. What remains is to calculate the $`a`$’s and $`\gamma `$’s in Eqs. (22)-(25). The results of these calculations are presented appendix B.

In principle it is possible to obtain the rate for all parameter values but to make the discussion more transparent we will consider regimes of the applies magnetic field $`(i)`$ $`E_BE_nN_{\text{QD}}^{1/2}`$ and $`(ii)`$ $`E_BE_nN_{\text{QD}}^{1/2}`$. In regime $`(i)`$ both $`\mathrm{\Delta }_0\mathrm{}c_\nu \mathrm{}^1`$ and $`\mathrm{\Delta }_0kT`$ (for experimentally relevant temperatures) and only the lowest order terms in $`\mathrm{\Delta }_0/\mathrm{}c_\nu \mathrm{}^1`$ need to be considered. In GaAs $`\mathrm{}c_l\mathrm{}^1=3.3\mathrm{}^1`$ nm$`\times `$meV and $`\mathrm{}c_t\mathrm{}^1=2.0\mathrm{}^1`$ nm$`\times `$meV for the longitudinal and transverse branches respectively. In the other regime we assume that the applied field dominates and that the internal field may be ignored, but no additional assumption are made in this case. This will in general require some numerical work. The resulting hyperfine-mediated spin-flip rates are

$`\mathrm{\Gamma }_{\text{HF}}`$ $`=`$ $`0.34{\displaystyle \frac{G_{\text{corr}}}{I^2}}{\displaystyle \frac{E_n^2}{N\text{QD}(\mathrm{}\mathrm{\Omega })^2}}{\displaystyle \frac{(eh_{14}\mathrm{})^2kT}{8\pi \rho c^5\mathrm{}^4}}\left(E_B^2+K_0^2+2E_BK_0\mathrm{cos}\theta \right)\text{ for }E_BE_nN_{\text{QD}}^{1/2}`$ (26)
$`\mathrm{\Gamma }_{\text{HF}}`$ $`=`$ $`{\displaystyle \frac{G_{\text{corr}}}{I^2}}{\displaystyle \frac{E_n^2(n(E_B)+1)}{N\text{QD}(\mathrm{}\mathrm{\Omega })^2}}{\displaystyle \frac{(eh_{14}\mathrm{})^2E_B^3}{8\pi \rho c^5\mathrm{}^4}}\left(C_0(E_B)+\left({\displaystyle \frac{E_B}{\mathrm{}\mathrm{\Omega }}}\right)^2C_2(E_B)\right)\text{ for }E_BE_nN_{\text{QD}}^{1/2}`$ (27)
Note that the rates have different dependence on the emitted energy $`\mathrm{\Delta }_0`$. In regime $`(i)`$ it is proportional to $`\mathrm{\Delta }_0^2`$ due to the Bose statistics, and in regime $`(ii)`$ it is proportional to $`\mathrm{\Delta }_0^3`$ for $`\mathrm{\Delta }_0\mathrm{}c_\nu \mathrm{}^1`$. In Eq. (27) we introduce the functions $`C_0`$ and $`C_2`$ which contain the details of the higher orbitals and the anisotropy integrals. For low fields $`\mathrm{\Delta }_0\mathrm{}c_\nu \mathrm{}^1`$ these functions are constant. The value for the spin flip rate in Eq. (26), for some typical value of $`K_010^2`$ T, is very low $`\mathrm{\Gamma }_{\text{HF}}10^6`$ s. This results in a lifetime of days, which will be extremely difficult to measure. We have plotted the general spin-flip rate in Eq. (27) for different confining energies and temperatures in Figs. 2 and 3, for both in-plane and perpendicular applied magnetic field. In the case of the higher confining energy there is very small difference between the in-plane and perpendicular direction of the external magnetic field. For the lower confining energy there is a substantial difference between the two directions of magnetic field. In this case the approximation $`\mathrm{\Delta }_0\mathrm{}c_\nu \mathrm{}^1`$ is no longer good and the the $`\mathrm{\Delta }_0^3`$ dependence of the rate is changed by the $`C`$-functions. The values of the rates are quite small, depending on the applied field, being in $`1`$ s<sup>-1</sup> for $`T=4`$ K at $`B0.5`$ T for a confining frequency of $`\mathrm{}\mathrm{\Omega }_0=2`$ meV.

## V Discussion

Generally speaking spin-flip rates require external dissipation. This is why at small Zeeman splittings they will contain small factor reflecting the vanishing phonon density of states. For a spin-orbit ratekhaetskii01:27 , the Kramer’s degeneracy results in this small factor being proportional to $`E_B^5`$. Presence of nuclear spins violate the Kramer’s theorem. So that the hyperfine rate discussed in the present paper is proportional to $`E_B^3`$ and will dominate at sufficiently low fields.

Comparing the hyperfine rate in Eq. (27) to spin-orbit related rateskhaetskii01:27  $`\mathrm{\Gamma }_{\text{SO}}`$ we obtain that the hyperfine rate dominates for fields below $`B0.3`$ T (assuming $`z_0=10`$ nm and $`\mathrm{}\mathrm{\Omega }_0=4`$ meV).

The role of the internal field produced by the nuclei is that the spin-flip rate does not vanish even for in the absence of external magnetic field. The show that the minimum rate is rather small, corresponding to a relaxation time of the order of days. We believe that the internal field will play an important role when the full dynamics of the electron spin in the presence of the nuclear system is considered.

Our model should also be applicable to other polar semiconductors which have non-zero nuclear spin, e.g. InAs where the $`g`$-factor is much larger.

## VI Aknowledgement

One author (SIE) would like to thank Daniela Pfannkuche, Alexander Chudnovskii and Daniel Loss for fruitful discussions. This work is a part of the research program of the “Stichting vor Fundementeel Onderzoek der Materie (FOM)”

## Appendix A Average over nuclear spin pair

When the thermal average over a pair of free nuclei the following nuclear correlation function appear in Eq. (21)

$`G_{\text{corr}}`$ $`=`$ $`(S_+^\alpha )^{}S_+^\beta \delta I_k^\alpha \delta I_k^{}^\beta _T`$ (28)
$`=`$ $`(S_+^\alpha )^{}S_+^\beta \delta \widehat{I}^\alpha \delta \widehat{I}^\beta _T\delta _{k,k^{}},`$
where $`\delta \widehat{I}^\alpha =\widehat{I}^\alpha I^\alpha _T`$ and the electron spin matrix elements being $`S_+^\alpha =𝒏_+|S^\alpha |𝒏_{}`$. The kronekker delta reflects that there are no correlations between two different nuclei and we have dropped the $`k`$ subscript in $`\delta \widehat{I}^\alpha \delta \widehat{I}^\beta _T`$ since the nuclei is assumed to be identical. By defining the symmetric correlator

$$g^{\alpha \beta }=\frac{1}{2}\delta \widehat{I}^\alpha \delta \widehat{I}^\beta +\delta \widehat{I}^\beta \delta \widehat{I}^\alpha _T$$
(29)
we get the following

$$G_{\text{corr}}=(S_+^\alpha )^{}S_+^\beta (g^{\alpha \beta }+i/2ϵ^{\alpha \beta \gamma }\widehat{I}^\gamma _T).$$
(30)
In an isotropic system, $`\widehat{𝑰}_T=0`$, the value of the correlation function is $`G_{\text{corr}}=\frac{1}{2}\frac{1}{3}I(I+1)`$.

## Appendix B Matrix elements for a parabolic quantum dot

Using the Darwin-Fock solutions and the Fang-Howard variational solution for the triangular quantum well we obtain the following equation are Eqs. (22)-(25)

$`a_{ll^{}}\gamma _{ll^{}}`$ $`=`$ $`\delta _{M,M^{}}{\displaystyle \frac{A^2Cn}{V_{\text{QD}}\zeta }}{\displaystyle \frac{(eh_{14}\mathrm{})^2\mathrm{\Delta }^3}{8\pi \rho c^5\mathrm{}^4}}(n(\mathrm{\Delta })+1){\displaystyle \frac{\mathrm{\Gamma }(n+n^{}+|M|+1)2^{3(n+n^{}+|M|)}}{n!n^{}!(n+|M|)!(n^{}+|M|)!}}`$ (31)
$`\times \left({\displaystyle \underset{\nu }{}}{\displaystyle \frac{c^5}{c_\nu ^5}}\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{}}{\mathrm{}c_\nu }}\right)^{2(n+n^{}+|M|1)}{\displaystyle _0^\pi }d(\mathrm{cos}\theta )A_\nu (\theta ){\displaystyle \frac{(\mathrm{sin}\theta )^{2(n+n^{}+|M|)}\mathrm{exp}\mathbf{\left(}\frac{1}{2}\mathbf{\left(}\frac{\mathrm{\Delta }\mathrm{sin}\theta }{\mathrm{}c_\nu \mathrm{}^1}\mathbf{\right)}^2\mathbf{\right)}}{\left(1+\left(\frac{\mathrm{\Delta }\mathrm{cos}\theta }{3\mathrm{}c_\nu z_0^1}\right)^2\right)^3}}\right)`$
where $`\zeta ^1=z_0𝑑z|\chi (z)|^4`$ and $`c^5=c_l^5+c_t^5`$ is the effective sound velocity of the phonons. The anisotropy functions are

$`A_t(\theta )`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}^2\theta (8\mathrm{cos}^4\theta +\mathrm{sin}^4\theta )}{4}}`$ (32)
$`A_l(\theta )`$ $`=`$ $`{\displaystyle \frac{9\mathrm{cos}^2\theta \mathrm{sin}^4\theta }{2}}.`$ (33)
The equation for $`\stackrel{~}{a}_{ll^{}}\stackrel{~}{\gamma }_{ll^{}}`$ is identical except for a different Kronekker delta function $`\delta _{M,M^{}}`$. The integral needs in general to be evaluated numerically but when $`\mathrm{\Delta }\mathrm{}c_\nu \mathrm{}^1`$ the exponential term and the denominator become unity and the resulting integral is simple to calculate.


text:
# Mel’nikov method revisited

## 1 Introduction

It is certainly impossible to give a fairly complete list of the papers devoted to the applications of the Mel’nikov method  for evaluating the onset of chaos arising in perturbed homoclinic (or heteroclinic) orbits. We will quote only some of the papers which are more directly connected with the present approach.

The main purpose of this paper is to illustrate a completely analytic procedure, based on a suitable extension of a classical method , which is parallel and – at least in part – complementary to the standard approach (see e.g. ). This procedure can be also applied to some “degenerate” situations, as to the case of unstable nonhyperbolic points, or of critical points located at the infinity (thus giving rise to unbounded orbits, e.g. the Keplerian parabolic orbits), and it is naturally “compatible” with the presence of general symmetry properties of the problem. For these reasons, apart from a clear “unifying” aspect, this method could be of great interest in celestial mechanics, and it could be a contribution to the study of some of the questions about the onset of chaos in unbounded phase space systems and in the presence of unbounded orbits. It has been remarked indeed that in this situation chaos manifests itself in a particularly dramatic way .

Let us remark immediately that, in the above mentioned degenerate cases, i.e. in the lack of the hypothesis of hyperbolicity, standard results of perturbation theory cannot be directly applied; for instance, to preserve the criticality, we will have to impose a sufficiently rapid vanishing of the perturbation at the critical point. We can then extend the introduction of Mel’nikov functions, and show not only the existence of smooth solutions of the perturbed problem, approaching the critical points and playing in this context the role of stable and unstable manifolds, but also the possible presence of infinitely many intersections of these asymptotic sets on the Poincaré sections, thus leading to a complicate dynamics typical of the homoclinic chaotic behaviour .

Being mainly interested in the methodological aspects, we will not devote special emphasis on new applications, but rather we will show, in the two last sections of this paper, how the method can be concretely applied in some typical situations, arising especially in celestial mechanics and general relativity.

## 2 Statement of the Problem

Although the procedure is completely general (indeed, we will state some results in a quite general setting in sections 3 and 4), we have actually in mind applications to celestial mechanics or to gravitational problems in general relativity; therefore we will restrict our attention mainly to problems described by Hamiltonians of the classical form

$$H=\frac{1}{2}\left(p_r^2+\frac{p_\theta ^2}{r^2}+\frac{p_\phi ^2}{r^2\mathrm{sin}^2\theta }\right)+V(r)$$
(1)
where $`V(r)`$ is some “long range” potential. Concretely, we have in mind potentials of the form

$$V(r)=\frac{1\beta _1}{r}\frac{\beta _2}{r^2}\frac{\beta _3}{r^3}$$
(2)
where $`\beta _i`$ are given parameters, which model several interesting situations, including standard Kepler problem, Kepler problems plus quadrupolar effects , Manev problem , or the motion of a charged particle in the field of a Schwarzschild black hole in general relativity , and many other situations.

Let us write the equations of the motion, assuming for a moment that it lies in the plane $`\theta =\pi /2`$,

$$\dot{r}=p_r,\dot{p}_r=\frac{dV}{dr}+\frac{L^2}{r^3},\dot{\phi }=\frac{L}{r^2},\dot{p}_\phi =0$$
(3)
where we have put $`L=p_\phi =`$const$`0`$. We are interested in this paper in the appearance of chaotic behaviour related to the presence of homoclinic (or heteroclinic) orbits subjected to perturbations; therefore, the relevant situations which may occur, depending on the values of the parameters $`\beta _i`$ in the potential (2), and which we are going to deal with in this paper, are the following:

i) the presence of an unstable equilibrium point $`r_u0`$ for the first two equations in the system (3), which involve the variables $`r,p_r`$. This point corresponds in the plane $`\theta =\pi /2`$ to a unstable circular orbit $`\gamma `$ of radius $`r_u`$. In this case, we have also the presence of a $`1`$parameter family of homoclinic bounded orbits biasymptotic to $`\gamma `$ (but see also case iii) below).

ii) the degenerate situation where the unstable equilibrium point is located at the infinity, i.e. $`r=\mathrm{},\dot{r}=0`$; the homoclinic orbits are in this case a family of parabolas.

iii) a “critical” case, with an unstable equilibrium point $`r_u0`$ and another unstable equilibrium point located at $`r=\mathrm{}`$, and a family of heteroclinic orbits connecting these points.

It is clear that – due to the spherical symmetry of the problem – all conclusions and properties stated for the plane $`\theta =\pi /2`$ are equally true for any plane for the origin in $`𝐑^3`$.

Let us now choose and single out the following homoclinic (heteroclinic in case iii)) orbit in the plane $`\theta =\pi /2`$, denoted by $`\widehat{\chi }(t)`$, written in the spherical variables

$$u:=(r,p_r,\phi ,p_\phi ,\theta ,p_\theta )$$

as follows

$$\widehat{\chi }(t):=(R(t),\dot{R}(t),\mathrm{\Phi }(t),L,\pi /2,0)$$
(4)
where $`R(t)`$ and $`\mathrm{\Phi }(t)`$ solve (3) with the conditions $`R(\pm \mathrm{})=r_u`$ in case i) or respectively $`R(\pm \mathrm{})=\mathrm{}`$ in case ii), and $`R(\mathrm{})=r_u,R(+\mathrm{})=+\mathrm{}`$ in case iii), and with $`\mathrm{\Phi }(0)=\pi `$. It is not necessary, for our purposes, to know explicitly the expression of the functions $`R(t)`$ and $`\mathrm{\Phi }(t)`$; it will be useful only to know that choosing $`R(0)=r_0`$ (the turning point) in cases i) and ii), $`R(t)`$ is an even function and $`\mathrm{\Phi }(t)`$ an odd function of the time $`t`$. Let us remark that any other homoclinic orbit can be transformed by means of a rotation into the $`\widehat{\chi }(t)`$ given by (4).

We now introduce a smooth (analytic) perturbation depending in general on all the variables $`u`$ and time-periodic; the equations of the motion we are considering are then, in general,

$`\dot{r}`$ $`=`$ $`p_r+ϵg_r(u,t)`$
$`\dot{p}_r`$ $`=`$ $`{\displaystyle \frac{p_\theta ^2}{r^3}}+{\displaystyle \frac{p_\phi ^2}{r^3\mathrm{sin}^2\theta }}{\displaystyle \frac{dV}{dr}}+ϵg_{p_r}(u,t)`$
$`\dot{\phi }`$ $`=`$ $`{\displaystyle \frac{p_\phi }{r^2\mathrm{sin}^2\theta }}+ϵg_\phi (u,t)\text{}`$ (5)
$`\dot{p}_\phi `$ $`=`$ $`ϵg_{p_\phi }(u,t)\text{}`$
$`\dot{\theta }`$ $`=`$ $`{\displaystyle \frac{p_\theta }{r^2}}+ϵg_\theta (u,t)`$
$`\dot{p}_\theta `$ $`=`$ $`{\displaystyle \frac{p_\phi ^2\mathrm{cos}\theta }{r^2\mathrm{sin}^3\theta }}+ϵg_{p_\theta }(u,t)`$
where $`ϵ1`$. Let us also write (5) in a more convenient compact form

$$\dot{u}=f(u)+ϵg(u,t)$$
(6)
where

$$f=J_uH$$
(7)
and $`J`$ is the standard symplectic matrix. Given any homoclinic orbit $`\chi (t)`$ of the unperturbed problem, one has that $`\chi _{t_0}:=\chi (tt_0)`$ satisfies, for all $`t_0𝐑`$,

$$\frac{d\chi _{t_0}}{dt}=f(\chi _{t_0})$$
(8)
In order to find conditions ensuring the occurrence of intersections of stable and unstable manifolds of the critical point for the perturbed problem, and hence the appearance of chaotic behaviour, we follow a (suitably extended) procedure which has been first used in this context (to the best of our knowledge) by Chow, Hale and Mallet-Paret  in a problem with $`1`$ degree of freedom. We have first to look for smooth solutions near the homoclinic orbit; we then put

$$u(t)=\chi (tt_0)+v(tt_0)$$
(9)
Substituting in (6), we find for $`v(t)`$ the equation (with the time shift $`tt+t_0`$)

$`\dot{v}`$ $`=`$ $`A(t)v+ϵg(\chi (t),t+t_0)+\mathrm{higher}\mathrm{order}\mathrm{terms}\mathrm{in}v\mathrm{and}ϵ`$
$`:=`$ $`A(t)v+G(t,t_0,ϵ)`$
having used the shorthand notation $`G(t,t_0,ϵ)`$, and where

$$A(t):=(_uf)(\chi (t))$$
(11)
This equation, or more often its first-order approximation

$$\dot{v}=A(t)v+ϵg(\chi (t),t+t_0)$$
(12)
is called the variational equation associated to (6) and to the orbit $`\chi (t)`$.

All the solutions $`v(t)`$ of (2) can be written, in implicit form (see )

$$v(t)=v_h(t)+\mathrm{\Psi }(t)_{t_1}^t\mathrm{\Psi }^1(s)G(s,t_0,ϵ)𝑑s$$
(13)
where $`t_1`$ is arbitrary, $`v_h(t)`$ is any solution of the homogeneous linear problem

$$\dot{v}_h=A(t)v_h$$
(14)
and $`\mathrm{\Psi }`$ is a fundamental matrix of solutions of (14). We now have to look for the solutions of the homogeneous equation (14) and to control their behaviour for $`t\pm \mathrm{}`$. As well known , this equation admits some solutions which remain bounded for all $`t𝐑`$ and other solutions which diverge for $`t\pm \mathrm{}`$. The asymptotic behaviour is exponential in the case of hyperbolic unstable points, is like some power $`|t|^\mu `$ in the case of unstable point at $`r=\mathrm{},\dot{r}=0`$, and more in general of nonhyperbolic critical points.

An obvious bounded solution is given by $`d\chi /dt`$. To find other solutions, the following results, which we will state for convenience in a quite general form, will be useful.

## 3 General Results. The Role of Symmetry

Let us start with the following definition (cf. )

###### Definition 1

A vector field

$$S:=\sigma (u)=\underset{i=1}{\overset{n}{}}\sigma _i(u)\frac{}{u_i}$$
(15)
is said to be a Lie-point symmetry (or – more exactly – the Lie generator of a Lie-point symmetry) for a general dynamical system (not necessarily Hamiltonian)

$$\dot{u}=f(u)$$
(16)
if

$$[\sigma ,f]=0$$
(17)
where $`[,]`$ denotes the usual Lie commutator.

We then can state:

###### Lemma 1

Let $`\chi (t)`$ be any given homoclinic orbit, solution of the problem (16), and assume that this problem admits a Lie-point symmetry $`S=\sigma (u)`$. Then

$$\zeta (t):=\sigma (\chi (t))$$

is a solution of the homogeneous part (14) of the variational equation associated to this problem.

Proof. From direct calculation, using (17) and (11),

$$\dot{\zeta }=\frac{d}{dt}\sigma (\chi (t))=_u\sigma \dot{\chi }=(f_u\sigma )(\chi (t))=(_uf\sigma )(\chi (t))=A(t)\zeta .$$

Remark. The obvious solution $`d\chi /dt`$ of (14) may be included into the solutions given by the above Lemma, indeed $`f=d/dt`$ is a (trivial) Lie-point symmetry of any autonomous dynamical system (16), expressing simply its time-flow invariance: if $`u_0(t)`$ is a solution, then $`u_0(t+t_0)`$ is also a solution.

###### Lemma 2

Given an Hamiltonian $`H`$, let $`K`$ be a constant of motion for the Hamiltonian, i.e. $`\{K,H\}=0`$. Then

$$S:=JK$$
(18)
is a Lie-point symmetry for the dynamical system $`\dot{u}=J_uH`$.

The proof is a straightforward verification.

## 4 General Results.
Mel’nikov-type Conditions

Still considering, from a general point of view, eq. (6) as arising from an Hamiltonian with $`n`$ degrees of freedom, assume that there are $`n`$ bounded solutions $`\zeta ^{[\alpha ]}(t)`$ and $`n`$ divergent solutions $`\eta ^{[\alpha ]}(t)`$ of the homogeneous part (14) of the variational equation. Let us construct the fundamental matrix $`\mathrm{\Psi }`$ in (13) putting the solutions $`\zeta ^{[\alpha ]}`$ in the first $`n`$ columns and the $`\eta ^{[\alpha ]}`$ in the remaining columns of $`\mathrm{\Psi }`$; recalling that $`det\mathrm{\Psi }=1`$ and observing that $`\mathrm{\Psi }^1`$ is a matrix having the $`\eta ^{[\alpha ]}`$ in the first rows and $`\zeta ^{[\alpha ]}`$ in the remaining rows, one easily sees (cf. ) that (13) can be written (still in implicit form, see (2))

$`v(t)=v_h(t)+`$ $`{\displaystyle \underset{\alpha }{}}`$ $`\zeta ^{[\alpha ]}(t){\displaystyle _{t_1}^t}\left(\eta ^{[\alpha ]}(s)JG(s,t_0,ϵ)\right)𝑑s`$
$`+`$ $`\eta ^{[\alpha ]}(t){\displaystyle _{t_1}^t}\left(\zeta ^{[\alpha ]}(s)JG(s,t_0,ϵ)\right)𝑑s`$
where $`v_h(t)`$ is a linear combination of the $`\zeta ^{[\alpha ]}(t),\eta ^{[\alpha ]}(t)`$.

We now have to look for the existence of bounded solutions of (2); more precisely we have to look for solutions $`v^{()}(t)`$ (and resp. $`v^{(+)}(t)`$) with the property of being bounded for $`t\mathrm{}`$ (resp. $`t+\mathrm{}`$); these will provide precisely those solutions

$$u^{(\pm )}(t)=\chi (tt_0)+v^{(\pm )}(tt_0)$$
(20)
of (6) which belong, by definition, to the unstable (resp. stable) manifold of the critical point.

If we linearize (4) around the solution $`v(t)0`$ (which amounts to deleting the higher-order terms in (2), or to considering the variational equation in the form (12)), and take into account the different asymptotic behaviour of the solutions $`\zeta ^{[\alpha ]}(t)`$ and $`\eta ^{[\alpha ]}(t)`$, it can be seen, as a consequence of the implicit-function theorem , that the existence of bounded solutions both for $`t\mathrm{}`$ and $`t+\mathrm{}`$ is ensured if the following Mel’nikov-type conditions are verified

$$M^{[\alpha ]}(t_0)=_{\mathrm{}}^+\mathrm{}\left(\zeta ^{[\alpha ]}(t)Jg(\chi (t),t+t_0)\right)𝑑t=0$$
(21)
Let us remark that if the perturbation is Hamiltonian, $`g=J_uW`$, as often happens in celestial mechanics, choosing $`\zeta ^{[1]}(t)=d\chi /dt`$, then the first of these conditions becomes

$$M^{[1]}(t_0)=_{\mathrm{}}^+\mathrm{}\{H,W\}(\chi (t),t+t_0)𝑑t=0$$
(22)
Similarly, if, e.g., $`\zeta ^{[2]}(t)`$ comes from a constant of motion $`K`$, according to Sect. 3, then

$$M^{[2]}(t_0)=_{\mathrm{}}^+\mathrm{}\{K,W\}(\chi (t),t+t_0)𝑑t=0$$
(23)
These conditions are identical to the conditions given e.g. in , where they are obtained by means of different procedures and hypotheses. It can be significant to remark that the present method then provides an extension of these formulas also to “degenerate” cases (nonhyperbolic points and possibly unbounded orbits), and to symmetries of more general nature, as in (15).

We have only to notice that, whereas the convergence of Mel’nikov integrals (21) is granted in the case of standard hyperbolic and isolated unstable points, the convergence is only “conditional”, i.e. along a suitable sequence of intervals (see  for any details), in cases i) and iii) of our classification in Sect. 2, and finally it must be controlled “by hand” in the non-hyperbolic case or in the case of critical point at the infinity. In the last cases one has to impose a sufficiently rapid vanishing of the contribution of the perturbation $`g(\chi (t),t)`$ when $`t\pm \mathrm{}`$, i.e., as expected, when approaching the critical point. The precise rate of this vanishing will depend on the specific problem in consideration (see Examples in the next sections, and also  for examples in $`1`$ degree of freedom). For what concerns the regularity of the solutions and of the asymptotic manifolds, see e.g. .

Changing now the point of view, and considering the Poincaré sections of the $`u^{()}`$ and $`u^{(+)}`$ solutions, the above arguments show that, once conditions (21) are satisfied, there occurs a crossing of the negatively and positively asymptotic sets on the Poincaré section . One usually imposes that the intersection is transversal; actually this condition is not strictly necessary, indeed it can be shown that it is sufficient that the crossing is “topological” , i.e., roughly, that there is really a “crossing”, from one side to the other, but we do not insist on this point, which goes beyond the scope of the present paper.

Thanks to the periodicity of the perturbation, one immediately deduces  that there is an infinite sequence of intersections of the positively and negatively asymptotic sets of the critical point in the Poincaré section, leading to a situation typical of the homoclinic chaos. The presence of such infinitely many intersections is clearly reminiscent of the chaotic behaviour expressed by the Birkhoff-Smale theorem in terms of the equivalence to the symbolic dynamics of Smale horseshoes. Actually, this theorem cannot be directly used in the present context because its standard proof is intrinsically based on hyperbolicity properties . However, several arguments can be invoked even in the “degenerate” cases, which allow us to conclude that, if the conditions (21) are satisfied, the perturbed problem exhibits a chaotic behaviour. We can refer e.g. to the classical arguments used in , and reconsidered by many others (see e.g. ), possibly resorting to singular coordinate transformations, such as the McGehee transformation in the case of critical point at the infinity, or the “blowing-up” method . More specifically, an equivalence to a “nonhyperbolic horseshoe” has been proved in , in which the contracting and expanding actions are not exponential but “polynomial” in time. The presence of Smale horseshoes and of a positive topological entropy has been also proved by means of a quite general geometrical or “topological” procedure  which holds, in the presence of area-preserving perturbations, even in the nonhyperbolic case.

## 5 Applications to Celestial Mechanics.
Two Degrees of Freedom

Coming back to the initial problems as stated in Sect. 2, the first step, according to the above discussion, is to write down the appropriate Mel’nikov conditions. It is clear that, apart from a rotation (this will change the expression of the perturbation $`g`$ in equations (5,6) into some new $`\stackrel{~}{g}`$ which will also depend in general on the parameters of the rotation: this will be explained in detail later), we can always assume that the homoclinic orbit lies in the plane $`z=0`$ and in particular is given precisely by the orbit $`\widehat{\chi }(t)`$ defined in (4), Sect. 2. This greatly simplifies the variational equation, and in particular its homogeneous part (14); indeed it is immediate to verify that, with this choice, the two equations for the variations $`v_5,v_6`$ of $`\theta `$ and $`p_\theta `$ are separated from the first four equations and admit regular bounded solutions; therefore the problem turns out to be 4-dimensional. An obvious symmetry for it is given by the rotations around the $`z`$axis, generated by

$$S=\frac{}{\phi }$$
(24)
Then, as discussed in Sect. 3, two bounded solutions of the homogeneous part of the variational equation are

$`\zeta ^{[1]}`$ $`=`$ $`{\displaystyle \frac{d\widehat{\chi }}{dt}}=(\dot{R}(t),\ddot{R}(t),\dot{\mathrm{\Phi }}(t),0,0,0)`$ (25)
$`\zeta ^{[2]}`$ $`=`$ $`(0,0,1,0,0,0)`$
where $`R(t),\mathrm{\Phi }(t)`$ have been defined in Sect. 2 (see (3,4)). We then get from (21) the two Mel’nikov conditions

$`M^{[1]}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}[\dot{R}(t)\stackrel{~}{g}_{p_r}(\widehat{\chi }(t),t+t_0)\ddot{R}(t)\stackrel{~}{g}_r(\mathrm{})`$
$`+\dot{\mathrm{\Phi }}(t)\stackrel{~}{g}_{p_\phi }(\mathrm{})]dt=0`$
$`M^{[2]}`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\stackrel{~}{g}_{p_\phi }(\widehat{\chi }(t),t+t_0)𝑑t=0`$ (27)
Notice that usually one has $`\stackrel{~}{g}_r=0`$. In the case where the problem, included the perturbation, is completely planar, and that the perturbation is generated by an Hamiltonian $`W`$ of the form

$$W=W(r,\phi ,t)$$
(28)
then, under rotation, $`W`$ is changed simply into $`W(r,\phi +\phi _0,t)`$ and the above conditions become

$$M^{[1]}(t_0,\phi _0)=_{\mathrm{}}^+\mathrm{}\left[\dot{R}(t)\frac{W(R(t),\mathrm{\Phi }(t)+\phi _0,t+t_0)}{r}+\dot{\mathrm{\Phi }}(t)\frac{W(\mathrm{})}{\phi }\right]𝑑t=0$$
(29)
and

$$M^{[2]}(t_0,\phi _0)=_{\mathrm{}}^+\mathrm{}\frac{W(R(t),\mathrm{\Phi }(t)+\phi _0,t+t_0)}{\phi }𝑑t=0$$
(30)
and it is easily seen that these can be also written according to the general form as given in (22,23).

As a first simple application of this situation, consider the classical Kepler problem with unstable equilibrium point at $`r=\mathrm{}`$ and $`\dot{r}=0`$ and a perturbation not depending on $`\phi `$, e.g.

$$W=\frac{\mathrm{sin}2\pi \nu t}{r^\delta }(\delta >1/2)$$
(31)
as in the classical Gyldén problem . Then, condition (30) is automatically satisfied, and the integral in (29) becomes (thanks to the parity of the function $`R(t)`$)

$$\mathrm{cos}2\pi \nu t_0_{\mathrm{}}^+\mathrm{}\frac{\dot{R}(t)}{R(t)^{\delta +1}}\mathrm{sin}2\pi \nu tdt$$

which clearly admits simple zeroes (the integral is easily seen to be $`0`$). In this example, it is also simple to evaluate explicitly the asymptotic behaviour of the solutions of the variational equation: one has indeed

$$\zeta =\dot{R}(t)|t|^{1/3},\eta (t)|t|^{4/3}\mathrm{for}|t|+\mathrm{}$$
(32)
and all conditions for the procedure given in Sect. 4 are satisfied, with – in particular – the condition $`\delta >1/2`$ which ensures in this case the correct rate of vanishing of the perturbation at $`r=\mathrm{}`$. Notice that this same condition would also guarantee that under the McGehee transformation , the perturbation is not singular. Then chaotic behaviour is expected for this problem.

A similar result holds for the more general (time-periodic) perturbations of the form (28) occurring e.g. in the restricted 3-body problems . These cases can be conveniently dealt with in this way. Assuming that $`W(R(t),\mathrm{\Phi }(t),t)0`$ for $`t\pm \mathrm{}`$, then (29) (or (22)) becomes

$$M^{[1]}(t_0,\phi _0)=_{\mathrm{}}^+\mathrm{}\frac{W}{t}(R(t),\mathrm{\Phi }(t)+\phi _0,t+t_0)𝑑t=0$$
(33)
(where clearly the derivative $`/t`$ must be performed only with respect to the explicit time-dependence of $`W`$); introducing then the “Mel’nikov potential” $`𝒲=𝒲(\phi _0,t_0)`$ (cf. ), corresponding to the perturbation $`W(r,\phi ,t)`$:

$$𝒲(t_0,\phi _0):=_{\mathrm{}}^+\mathrm{}W(R(t),\mathrm{\Phi }(t)+\phi _0,t+t_0)𝑑t$$
(34)
then one gets from this definition and from (33,30)

$$M^{[1]}(t_0,\phi _0)=\frac{𝒲}{t_0}=0,M^{[2]}(t_0,\phi _0)=\frac{𝒲}{\phi _0}=0$$
(35)
In other words, the two Mel’nikov conditions are equivalent to the existence of stationary points for the Mel’nikov potential $`𝒲(\phi _0,t_0)`$. On the other hand, $`𝒲`$ is a smooth doubly-periodic function, and such a function certainly possesses points $`\overline{t_0},\overline{\phi _0}`$ where the two partial derivatives in (35) vanish, and this implies that the two conditions (33,30) are certainly satisfied (one has only to check that these stationary points of $`𝒲`$ are isolated).

## 6 Applications to Celestial Mechanics.
Three Degrees of Freedom

Let us consider finally a perturbation depending on both angles $`\phi ,\theta `$. As already stated, to obtain Mel’nikov condition in the form (5,27), one has to transform the generic homoclinic orbit into the orbit $`\widehat{\chi }(t)`$ given by (4). This is obtained by means of a rotation defined by the following Euler angles (with the conventions and notations as in ):

$$\mathrm{\Omega },i,\omega $$
(36)
where, with the language of celestial mechanics, $`i`$ is the inclination of the plane of the orbit, $`\omega `$ the angle of the perihelion with the line of nodes in the orbital plane, and $`\mathrm{\Omega }`$ is the longitude of the ascending node.

Considering for simplicity a perturbation generated by an Hamiltonian not depending on the variables $`p`$, which we denote here, with a little abuse of notation, either by $`W(r,\phi ,\theta ,t)`$ or by $`W(𝐱,t),𝐱=(x,y,z)`$, and denoting by $`B`$ the matrix of this rotation, the new expression $`\stackrel{~}{W}`$ of the perturbation is obtained by replacing $`𝐱`$ with $`B𝐱`$. It is then easy to verify that Mel’nikov conditions (5,27) become

$`M^{[1]}(t_0,\omega ,\mathrm{\Omega })`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left[\dot{R}{\displaystyle \frac{W}{r}}+\dot{\mathrm{\Phi }}\left({\displaystyle \frac{\stackrel{~}{W}}{\phi }}\right)_0\right]𝑑t=0`$ (37)
$`M^{[2]}(t_0,\omega ,\mathrm{\Omega })`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}\left({\displaystyle \frac{\stackrel{~}{W}}{\phi }}\right)_0𝑑t=0`$ (38)
where

$$\left(\frac{\stackrel{~}{W}}{\phi }\right)_0=R(t)\left(C_1(R(t),\mathrm{\Phi }(t),\omega ,\mathrm{\Omega })\mathrm{sin}\mathrm{\Phi }(t)+C_2(\mathrm{})\mathrm{cos}\mathrm{\Phi }(t)\right)$$
(39)
with

$`C_1`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{W}}{x_1}}\right)_0(\mathrm{cos}\omega \mathrm{cos}\mathrm{\Omega }\mathrm{cos}i\mathrm{sin}\omega \mathrm{sin}\mathrm{\Omega })+\text{}`$
$`+`$ $`\left({\displaystyle \frac{\stackrel{~}{W}}{x_2}}\right)_0(\mathrm{cos}\omega \mathrm{sin}\mathrm{\Omega }\mathrm{cos}i\mathrm{sin}\omega \mathrm{cos}\mathrm{\Omega })+\left({\displaystyle \frac{\stackrel{~}{W}}{x_3}}\right)_0\mathrm{sin}i\mathrm{sin}\omega `$
$`C_2`$ $`=`$ $`\left({\displaystyle \frac{\stackrel{~}{W}}{x_1}}\right)_0(\mathrm{sin}\omega \mathrm{cos}\mathrm{\Omega }+\mathrm{cos}i\mathrm{cos}\omega \mathrm{sin}\mathrm{\Omega })+`$
$`+`$ $`\left({\displaystyle \frac{\stackrel{~}{W}}{x_2}}\right)_0(\mathrm{sin}\omega \mathrm{sin}\mathrm{\Omega }+\mathrm{cos}i\mathrm{cos}\omega \mathrm{cos}\mathrm{\Omega })\left({\displaystyle \frac{\stackrel{~}{W}}{x_3}}\right)_0\mathrm{sin}i\mathrm{cos}\omega `$
and where $`(\stackrel{~}{W}/x_i)_0`$ means that in the derivative of the given $`W`$ with respect to $`x_i`$ one has to replace $`𝐱`$ with $`B𝐱`$ and finally put $`z=0`$ (or $`\theta =\pi /2)`$.

If $`i=0`$, i.e. if the problem is completely planar, including the perturbation, or if the perturbation is “generic”, i.e. has no “preferred” direction in the space (as often happens, see below for an example), and therefore it is not restrictive to choose $`i=0`$, then one gets

$$C_1=\left(\frac{\stackrel{~}{W}}{x_1}\right)_0\mathrm{cos}\phi _0+\left(\frac{\stackrel{~}{W}}{x_2}\right)_0\mathrm{sin}\phi _0,C_2=\left(\frac{\stackrel{~}{W}}{x_1}\right)_0\mathrm{sin}\phi _0\left(\frac{\stackrel{~}{W}}{x_2}\right)_0\mathrm{cos}\phi _0$$

where $`\phi _0=(\omega +\mathrm{\Omega })`$, and the above expressions (29,30) are recovered.

An example is provided by a problem in general relativity . Consider indeed the motion of a relativistic charged particle in a gravitational field produced by a Schwarzschild black hole, perturbed by a homogeneous constant electric field. It can be shown  that the perturbation is given by

$$W=(r)(l_1\mathrm{cos}\phi \mathrm{sin}\theta +l_2\mathrm{sin}\phi \mathrm{sin}\theta +l_3\mathrm{cos}\theta )$$
(40)
where $`(r)`$ is suitable function, and $`𝐥=(l_1,l_2,l_3)`$ is the direction of the perturbing field. This direction is generic and therefore no rotation is required (it would simply change the values of $`l_i`$, which are not fixed).

In this case, the first Mel’nikov condition in (37) is identically satisfied, because the perturbation is independent of the time (cf. (33)), and the other one becomes

$$M^{[2]}=M^{[2]}(\phi _0)=(l_1\mathrm{sin}\phi _0l_2\mathrm{cos}\phi _0)_{\mathrm{}}^+\mathrm{}(R(t))\mathrm{cos}\mathrm{\Phi }(t)𝑑t$$
(41)
which clearly admits simple zeroes, if only one can show that the integral is different from zero (see ).

Notice that in this example one has that the effective unperturbed potential admits bounded homoclinic orbits biasymptotic to an unstable point $`r_u0`$ (which corresponds to unstable circles, as in case i) of our initial classification), therefore the above integral is expected to converge only conditionally .

Similar results hold if the perturbation is produced by a constant homogeneous magnetic field . It can be observed that, whereas the component of the electric and of the magnetic field on the plane of motion leads to a chaotic dynamics, the component normal to the plane does not. Since the problem is spherically symmetric, this argument can be applied to every plane for the origin. Thus, given an electric or magnetic field, on each plane for the origin (except at most the one normal to the field), chaos appears for a suitable choice of initial conditions.


text:
# Pairing of Cooper Pairs in a Fully Frustrated Josephson Junction Chain

\[

## Abstract

We study a one-dimensional Josephson junction chain embedded in a magnetic field. We show that when the magnetic flux per elementary loop equals half the superconducting flux quantum $`\varphi _0=h/2e`$, a local $`\mathrm{𝖹𝖹}_2`$ symmetry arises. This symmetry is responsible for a nematic Luttinger liquid state associated to bound states of Cooper pairs. We analyze the phase diagram and we discuss some experimental possibilities to observe this exotic phase.

\]

During the last twenty years, Josephson junction arrays have proved to be very good tools to investigate classical and quantum phase transitions . Recently, much attention has been payed to systems which display highly degenerate classical ground states  due to the presence of Aharonov-Bohm cages. Interestingly, a glassy vortex phase without disorder has been predicted for such two-dimensional ($`2D`$) structures  in agreement with experimental observations . In the present work, we investigate the influence of quantum fluctuations on such systems in a $`1D`$ model introduced in . In this remarkably simple example, the huge classical degeneracy is a direct consequence of a local $`\mathrm{𝖹𝖹}_2`$ symmetry which is unbroken in the presence of quantum fluctuations. We show that this may stabilize an unusual nematic Luttiger liquid (LL) phase in which charge $`4e`$ bound states of Cooper pairs are the elementary objects.

We consider the chain of loops shown in Fig. 1 embedded in a uniform magnetic field. We denote by $`\varphi `$ the magnetic flux per elementary plaquette and we set $`\gamma =2\pi \varphi /\varphi _0`$ where $`\varphi _0=h/2e`$ is the superconducting flux quantum. Each site of this lattice is supposed to be occupied by a superconducting island. A convenient description of the low-energy Hilbert space of this system involves local boson operators $`a_n^{},b_n^{},c_n^{}`$ $`(a_n,b_n,c_n)`$ that create (destroy) Cooper pairs on the three types of islands of the lattice, respectively represented by black, grey and white circles in Fig. 1. The system is described by the following Josephson coupling Hamiltonian :

$$H_\mathrm{J}=t_\mathrm{J}\underset{n}{}a_n^{}(b_n+c_n+b_{n1}+\mathrm{e}^{\mathrm{i}\gamma }c_{n1})+\text{h.\hspace{0.17em}c.}$$
(1)
We first focus on the special value $`\gamma =\pi `$ (half a flux quantum per loop). As shown in , this Hamiltonian has, in this case, a single particle spectrum composed of three highly degenerate flat bands $`\epsilon _0=0`$, $`\epsilon _\pm =\pm 2t_\mathrm{J}`$. The corresponding eigenstates can be chosen as strictly localized (cage states) around each fourfold coordinated site (see Fig. 1). This leads naturally to the notion of Aharonov-Bohm cages discussed in .

Let us introduce the set of boson operators $`A_{\alpha ,n}^{}`$ ($`A_{\alpha ,n}`$) that create (destroy) one Cooper pair with energy $`\epsilon _\alpha `$ ($`\alpha =0,\pm `$) in a cage state localized around the $`n`$th fourfold site. These operators can be simply expressed as a linear combination of the operators $`a_n^{},b_n^{},c_n^{},b_{n1}^{},c_{n1}^{}`$ only, whose coefficients are given in Fig. 1 so that we get :

$$H_\mathrm{J}=\underset{n,\alpha }{}\epsilon _\alpha A_{\alpha ,n}^{}A_{\alpha ,n}\text{ . }$$
(2)
In this present form, $`H_\mathrm{J}`$ clearly exhibits a local $`U(1)`$ symmetry. We shall now study the effect of boson-boson interaction on this symmetry. Therefore, we consider a real valued function $`ns_n`$, and we construct a unitary operator $`U_s`$ defined by : $`U_sA_{\alpha ,n}U_s^1=\mathrm{e}^{\mathrm{i}s_n}A_{\alpha ,n}`$ which commutes with $`H_\mathrm{J}`$. Using the precise form of the cage states, we easily obtain :

$`U_sa_n^{}a_nU_s^1`$ $`=`$ $`a_n^{}a_n`$ (3)
$`U_sb_n^{}b_nU_s^1`$ $`=`$ $`\mathrm{cos}^2(\mathrm{\Delta }_n)b_n^{}b_n+\mathrm{sin}^2(\mathrm{\Delta }_n)c_n^{}c_n+z_n`$ (4)
$`U_sc_n^{}c_nU_s^1`$ $`=`$ $`\mathrm{sin}^2(\mathrm{\Delta }_n)b_n^{}b_n+\mathrm{cos}^2(\mathrm{\Delta }_n)c_n^{}c_nz_n\text{ ,}`$ (5)
where we have set $`z_n=\mathrm{i}\mathrm{sin}(2\mathrm{\Delta }_n)(b_n^{}c_nc_n^{}b_n)/2`$ and $`\mathrm{\Delta }_n=(s_{n+1}s_n)/2`$. From these transformation laws, we can readily see that any interaction term involving the bilinear operators $`a_m^{}a_m`$, $`b_n^{}b_n+c_n^{}c_n`$ preserves the local $`U(1)`$ symmetry. Physically, this symmetry implies that the total number of bosons in each cage is separately conserved and the system remains an insulator. However, this symmetry is fragile since it is easy to find two-body interactions which break it. For instance, a Hubbard-like interaction term $`_n(a_n^{}a_n)^2+(b_n^{}b_n)^2+(c_n^{}c_n)^2`$, has this effect which is manifested by the appearance of delocalized two-particle bound states discussed in . An exciting feature of this system is that this type of interaction still preserves a subgroup of the full $`U(1)`$ corresponding to a local $`\mathrm{𝖹𝖹}_2`$ symmetry. This subgroup corresponds to $`s_n=0[\pi ]`$ for all $`n`$. With this restriction, it is easy to check that the operator $`b_m^{}b_mb_n^{}b_n+c_m^{}c_mc_n^{}c_n`$ commutes with $`U_s`$ for all $`(m,n)`$. This local $`\mathrm{𝖹𝖹}_2`$ symmetry has an important physical consequence since it means that the parity of the total number of bosons in each cage is separately conserved. Therefore, if two-particle interactions lead to coherent transport through the chain for a many-boson system, quasi off-diagonal long range order may occur only for composite objects built with an even number of original bosons i. e., here, of Cooper pairs. In other words, a superconducting Josephson junction chain with this geometry and half a flux quantum per loop may realize a quasi-Bose condensate (in fact a LL) of charge $`4e`$ composite bosons !

To discuss in more details the physics of this system, it is useful to rephrase these symmetry considerations in the language of quantum rotor mo dels. These offer the advantage of an intuitively simple classical limit defined from the phase of superconducting order parameter. Formally, we introduce three phase fields $`\theta _n,\phi _n,\chi _n`$ and their canonically conjugate fields $`\mathrm{\Pi }_{\theta ,n},\mathrm{\Pi }_{\phi ,n},\mathrm{\Pi }_{\chi ,n}`$, which are related to the local Bose operators by : $`a_n^{}=\mathrm{\Pi }_{\theta ,n}^{1/2}\mathrm{e}^{\mathrm{i}\theta _n}`$, $`b_n^{}=\mathrm{\Pi }_{\phi ,n}^{1/2}\mathrm{e}^{\mathrm{i}\phi _n}`$, $`c_n^{}=\mathrm{\Pi }_{\chi ,n}^{1/2}\mathrm{e}^{\mathrm{i}\chi _n}`$. Assuming that the local particle number fluctuations are small, we get the quantum phase Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \frac{E_\mathrm{C}}{2}}{\displaystyle \underset{n}{}}\mathrm{\Pi }_{\theta ,n}^2+\mathrm{\Pi }_{\phi ,n}^2+\mathrm{\Pi }_{\chi ,n}^2`$ (8)
$`E_\mathrm{J}{\displaystyle \underset{n}{}}\mathrm{cos}(\theta _n\phi _n)+\mathrm{cos}(\theta _n\chi _n)+`$
$`\mathrm{cos}(\theta _n\phi _{n1})+\mathrm{cos}(\theta _n\chi _{n1}\gamma )\text{ ,}`$
where $`E_\mathrm{C}`$ is the charging energy and $`E_\mathrm{J}`$ the Josephson coupling between islands. Note that the present modelling of capacitive effects is not meant to be very realistic, since for the sake of simplicity, we have not taken into account off-diagonal elements of the capacitance matrix. This choice corresponds to a local Hubbard-like interaction term between Cooper pairs. For convenience, we set $`\sqrt{E_\mathrm{C}E_\mathrm{J}}=1`$.

The classical ground state of $`H`$ is easily obtained for any $`\gamma `$. Indeed, if we set $`x_n=\theta _{n+1}\theta _n\gamma /2`$, and eliminate $`\phi _n`$ and $`\chi _n`$, minimizing $`H`$ is equivalent to minimize $`F(x_n)=\left|\mathrm{cos}\left(x_n/2+\gamma /4\right)\right|\left|\mathrm{cos}\left(x_n/2\gamma /4\right)\right|`$ for all $`n`$. As shown in Fig. 2, $`F`$ has two local minima in $`x_n=0`$ and $`x_n=\pi `$. For $`0<\gamma <\pi `$, one has $`F(0)<F(\pi )`$ and the classical ground state is unique (up to a global $`U(1)`$ degeneracy). By contrast, for $`\gamma =\pi `$, one has $`F(0)=F(\pi )`$ so that, for a given plaquette, we obtain two degenerate ground states (up to a global translation of the phase variables) which are illustrated in Fig. 3.

These states only differ in the sign of the superconducting currents which circulate around the plaquette. For a chain made up of $`N`$ loops, we thus get $`2^N`$ degenerate classical ground states up to a global translation of the phase variables.

This huge degeneracy is a direct consequence of a local $`\mathrm{𝖹𝖹}_2`$ symmetry of $`H`$. Note that $`H`$ is not invariant under the full local $`U(1)`$ group related to the Aharonov-Bohm cages. This occurs since the Josephson term in $`H`$ may be written as a strongly non linear expression of the basic local Bose operators. For the $`\mathrm{𝖹𝖹}_2`$ transformations, it is an easy task to translate the $`U_s`$ operators in the language of phase variables. An interesting local $`\mathrm{𝖹𝖹}_2`$ transformation is provided by a kink in the $`s_m`$’s. Let $`U_n`$ be the transformation defined by $`s_m=0`$ for $`mn`$ and $`s_m=\pi `$ for $`m>n`$. This transformation does not modify the phase variables for $`mn`$ whereas it shifts them by $`\pi `$ if $`m>n`$ and it permutes $`\phi _n`$ and $`\chi _n`$. Thus, its main physical effect is to change the currents flowing around the plaquette located between $`n`$ and $`n+1`$ into their opposite value. From this description, we deduce that starting from a given classical ground state, we may generate any other ground state by applying a finite sequence of such $`U_n`$ operators. We also see that, in the classical limit considered here, the local $`\mathrm{𝖹𝖹}_2`$ symmetry is spontaneously broken, yielding ground states with well-defined local circulating supercurrents. Note that $`U_n`$ also leaves most conjugate variables unchanged except $`\mathrm{\Pi }_{\phi ,n}`$ and $`\mathrm{\Pi }_{\chi ,n}`$ which are exchanged. As a result, we may add to $`H`$ any term involving these conjugate fields without breaking the local $`\mathrm{𝖹𝖹}_2`$ symmetry, provided the spatial symmetry between the $`\phi _n`$ and $`\chi _n`$ degrees of freedom is respected. Experimentally, this would require a tight control of offset charges since these may seriously alter an otherwise excellent geometrical symmetry of the chain.

For real systems, it may become important to take into account quantum fluctuations of the phase variables, especially when the superconducting islands are so small that their charging energy $`E_\mathrm{C}`$ can no longer be neglected in comparison to the Josephson coupling energy $`E_\mathrm{J}`$. For a single loop and $`\gamma =\pi `$, these quantum fluctuations have been shown, theoretically  and experimentally , to induce tunneling between the two degenerate classical ground states shown in Fig. 3. The true quantum mechanical ground state is therefore a macroscopic linear superposition of these two classical states and provides a simple example of a “Schrödinger cat”. For a system with $`N`$ loops, eigenstates are classified according to the various irreducible representations of the local $`\mathrm{𝖹𝖹}_2`$ group which mix all the $`2^N`$ classical ground states. One of our next goals is to describe how quantum fluctuations lift the degeneracy among these representations, which is an artifact of the classical limit.

At small $`g=\sqrt{E_\mathrm{C}/E_\mathrm{J}}`$, the properties of the system are actually very similar to those of a quantum $`XY`$ model. For small $`\gamma `$, we thus expect the infinite chain to be in a LL phase for $`g<g^{}(\gamma )`$ and in a gapful insulating (I) phase for $`g>g^{}(\gamma )`$. The transition at $`g^{}(\gamma )`$ is of Berezinskii-Kosterlitz-Thouless (BKT) type . Simple spin wave calculations using the harmonic approximation of $`H`$ around its classical ground state predict $`g^{}(\gamma )=\sqrt{\mathrm{cos}(\gamma /4)}g^{}(0)`$ with $`g^{}(0)=\pi \sqrt{3}/2`$. The main effect of the magnetic field, in this simple approximation, is thus to replace $`g`$ by an effective $`g_{eff}=g/\sqrt{\mathrm{cos}(\gamma /4)}`$ which controls all the correlation function exponents in the LL phase.

To analyze the effect of the quantum fluctuations in the vicinity of $`\gamma =\pi `$ where the additional local $`\mathrm{𝖹𝖹}_2`$ symmetry emerges, it is convenient to eliminate the twofold coordinated islands to get a simple description of the low-energy physics of this system. Therefore, instead of $`H`$, we now consider the following Hamiltonian :

$`H_{XY}`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{g^{}}{2}}\mathrm{\Pi }_{\theta ,n}^2{\displaystyle \frac{1}{g^{}}}\{\mathrm{cos}[p(\theta _n\theta _{n+1}\gamma /2)]+`$ (10)
$`ϵ\mathrm{cos}(\theta _n\theta _{n+1}\gamma /2)\}\text{ .}`$
The parameter $`g^{}`$ is provided by fitting the exponent of the correlation function $`\mathrm{e}^{\mathrm{i}(\theta _m\theta _n)}`$ in the semi-classical regime with its value obtained with $`H`$ in the harmonic approximation. This choice leads to $`g^{}=2^{5/4}3^{1/2}g`$. The parameter $`ϵ=4|\gamma \pi |`$ is determined from the energy splitting between the two local minima of the single loop potential energy. Finally, in our case, we have $`p=2`$ but we discuss thereafter the properties of $`H_{XY}`$ for an arbitrary $`p`$.

The Hamiltonian $`H_{XY}`$ has a local $`\mathrm{𝖹𝖹}_p`$ symmetry at $`ϵ=0`$, corresponding to the local transformations $`T_a`$ : $`\theta _j\theta _j+2\pi a_j/p`$ where $`a_j`$ is an integer. The irreducible representations of this group are easily obtained in a basis which diagonalizes simultaneously the $`\mathrm{\Pi }_{\theta ,j}`$’s. For a state $`|\psi `$ such that $`\mathrm{\Pi }_{\theta ,j}|\psi =l_j|\psi `$, where $`l_j`$ is an integer, we have $`T_a|\psi =\mathrm{exp}(\mathrm{i}_j\frac{2\pi }{p}l_ja_j)|\psi `$. Writing $`l_j=m_j+pn_j`$ with $`m_j`$ and $`n_j`$ integers and $`0m_jp1`$, we find that the set of $`m_j`$’s completely specifies the irreducible representation of the local $`\mathrm{𝖹𝖹}_p`$ group. For each such representation, the action of the corresponding projector on the approximate gaussian ground state of $`H_{XY}`$ produces a natural trial wave-function at least when $`g^{}1`$. We have computed the expectation value of $`H_{XY}`$ on these states. Doing so, we noticed that $`2\pi /p`$-tunnel processes occuring on different lattice sites are mostly uncorrelated. Neglecting completely these correlations we get :

$$H_{XY}=\frac{g^{}L^2}{2N}\frac{C\mathrm{e}^f}{g^{}}\underset{j}{}\mathrm{cos}\left[\frac{2\pi }{p}(m_jL/N)\right]\text{,}$$
(11)
up to a constant energy independent of the representation and to factors of order $`\mathrm{e}^{2f}`$. In (11), $`L=_jl_j`$ is the total angular momentum, $`C`$ is a number close to $`2\pi ^28`$ at small $`g^{}`$, and $`f4\pi /pg^{}`$. The ground state is therefore obtained by choosing the identity representation of the local $`\mathrm{𝖹𝖹}_p`$ group ($`m_j=0`$).

Next, we see that the term proportionnal to $`ϵ`$ couples different irreducible representations of the local $`\mathrm{𝖹𝖹}_p`$ group. When $`p=2`$ the action of this perturbation on the $`2^N`$ low-energy trial states just discussed, is well described by a quantum Ising model in a transverse magnetic field :

$`H_\mathrm{I}`$ $`=`$ $`{\displaystyle \frac{1}{g^{}}}\left(C\mathrm{e}^f{\displaystyle \underset{n}{}}\sigma _n^X+Dϵ{\displaystyle \underset{n}{}}\sigma _n^Z\sigma _{n+1}^Z\right)\text{ ,}`$ (12)
where $`D`$ is close to $`1`$ for small $`g^{}`$. In terms of these Ising variables, the local $`\mathrm{𝖹𝖹}_2`$ symmetry corresponds to interchanging the $`|+`$ and $`|`$ states on any given subset of sites. It is therefore implemented by the $`\sigma _n^X`$ operators. This model has a continuous phase transition (in the universality class of the $`2D`$ Ising model) at its self-dual point which corresponds to $`ϵ=(C/D)\mathrm{e}^{2\pi /g^{}}`$. Furthermore, it is easy to see that $`H_{XY}`$ exhibits a BKT transition at $`\gamma =\pi `$ (so $`ϵ=0`$) corresponding to the loss of quasi long range order for the nematic order parameter $`\mathrm{e}^{2\mathrm{i}\theta }`$. In the harmonic approximation, this occurs at $`\stackrel{~}{g}(\pi )=g^{}(\pi )/4`$.

Gathering these informations obtained in various limits, we get the phase diagram drawn in Fig. 4. Besides the two familiar phases, namely the (I) phase at large $`g`$ and the LL phase characterized by an algebraic order of the $`\mathrm{e}^{\mathrm{i}\theta }`$ order parameter, the most interesting result is the presence of a remarkable nematic Luttinger liquid (NLL) phase which may be viewed as a quasi-ordered condensate of pairs of Cooper pairs associated to the order parameter $`\mathrm{e}^{2\mathrm{i}\theta }`$.

The dotted line of $`2D`$ Ising-type separates the two Luttinger phases. The Ising order parameter is vanishing in the NLL phase and builds up in the conventional LL phase where its square is proportional to the algebraically decaying part of the $`\mathrm{e}^{\mathrm{i}\theta }`$ autocorrelation function. The physical picture presented here has already been uncovered by Lee and Grinstein  in the framework of a $`2D`$ classical $`XY`$ model with squared cosine interaction. In their language, the BKT transition from the $`2e`$ LL phase to the I phase involves $`2\pi `$ vortices of the $`\theta `$ field, by contrast to $`\pi `$ vortices between the NLL phase and the I phase. The nature of the multicritical point $`Q`$ remains mysterious to us. Note that the coupling between $`XY`$ and Ising-like degrees of freedom has been extensively studied , already in the context of fully frustrated Josephson juntion arrays. Nevertheless, in these works, the Ising order parameter has a very different nature since it describes the possible melting of the vortex lattice. It implies that the $`XY`$ transition has to appear at lower temperatures than the Ising transition . In our case and in the vicinity of $`\gamma =\pi `$, the $`XY`$ transition occurs for larger $`g`$ than the Ising one. So, the effect of the Ising domain walls on the phase stiffness is very different in the two situations.

To conclude, we have established the possibility of macroscopic condensation (in the sense of a LL) of charge $`4e`$ objects in a Josephson junction chain for $`\gamma \pi `$. Experimentally, it may be possible to detect this binding of Cooper pairs by connecting such a network to superconducting leads. Indeed, we expect here a phenomenon analogous to Andreev’s reflection : an ordinary Cooper pair of charge $`2e`$ entering the chain at low energy (compared to the gap between different representations of the local $`\mathrm{𝖹𝖹}_2`$ group ) will leave behind a pair of charge $`2e`$ so that a charge $`4e`$ composite object may propagate along the chain. Another possibility is to close the chain into a large ring. In this geometry, we expect quantum oscillations of the global current with respect to the magnetic flux across the ring with an elementary period $`\varphi _0/2`$ as long as $`\gamma \pi `$.

We thank P. Lecheminant, D. Mouhanna, B. Pannetier and E. Serret for fruitful and stimulating discussions.