Characterization of the quantum phase transition in a two-mode Dicke model for different cooperation numbers
L. F. Quezada E. Nahmad-Achar







luis.fernando@correo.nucleares.unam.mx

Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de
México, Apartado Postal 70-543, 04510 Ciudad de México, México.




nahmad@nucleares.unam.mx

Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de
México, Apartado Postal 70-543, 04510 Ciudad de México, México.


We show how the use of variational states to approximate the ground
state of a system can be employed to study a multi-mode Dicke model.
One of the main contributions of this work is the introduction of
a not very commonly used quantity, the cooperation number, and the
study of its influence on the behavior of the system, paying particular
attention to the quantum phase transitions and the accuracy of the
used approximations. We also show how these phase transitions affect
the dependence of the expectation values of some of the observables
relevant to the system and the entropy of entanglement with respect
to the energy difference between atomic states and the coupling strength
between matter and radiation, thus characterizing the transitions
in different ways.

Characterization of the quantum phase transition in a two-mode Dicke
model for different cooperation numbers
    E. Nahmad-Achar
    Received: date / Accepted: date
============================================================================================================





§ INTRODUCTION


Quantum phase transitions (QPTs) are informally seen as sudden, drastic
changes in the physical properties of the ground state of a system
at zero temperature due to the variation of some parameter involved
in the modeling Hamiltonian. One model of particular interest for
the study of such phenomena is the Dicke model <cit.>, as it
describes, in a simplified way (electric dipole approximation), the
interaction between matter and electromagnetic radiation. In 1973,
Hepp and Lieb <cit.>, and Wang and Hioe <cit.>
first theoretically proved the existence of a second-order QPT in
the Dicke model. Wang and Hioe also treated the multi-mode radiation case, where they reduce it to a single-mode case by using an effective coupling constant. To date, this QPT has been experimentally observed
in a Bose-Einstein Condensate coupled to an optical cavity <cit.>
and it has been shown to be relevant to quantum information and quantum
computing <cit.>.

Even though the formal definition of a QPT requires us to compute
the ground state's energy as a function of any desired parameter in
order to find its transition values, one of the main contributions
of this work is to show how the QPT in the Dicke model influences
the behavior of other quantities relevant to the system, thus characterizing
the transition in different, simpler ways.




  §.§.§ Quantum Phase Transitions


The formal definition of the concept of “quantum phase” that we
will be using throughout this paper is that of an open region ℛ⊆ℝ^ℓ
where the ground state's energy ℰ_0, as a function
of ℓ parameters involved in the modeling Hamiltonian, is analytic.
Thus a QPT is identified by the boundary ∂ℛ of
the region at which ∂^nℰ_0/∂ x^n
is discontinuous for some n (known as the order of the transition).

Notice that in the previous definition, for the sake of generality,
we did not consider the thermodynamic limit, as it has been shown
that interesting phenomena regarding QPTs occur even for a finite
number of particles <cit.>.










  §.§.§ Modeling Hamiltonian


The Hamiltonian (Dicke's Hamiltonian) describing the interaction,
in a dipolar approximation, between N two-level identical atoms (same
energy difference between the two levels) and one-mode of an electromagnetic
field in an ideal cavity, has the expression (ħ=1)


    H_D=ω_AJ_z+Ω a^†a-γ/√(N)(J_-+J_+)(a+a^†).



Here, ω_A is the energy difference between the atomic levels,
Ω is the frequency of the field's mode, γ is the dipolar
coupling constant, J_z, J_-, J_+ are the collective
spin operators and a, a^† are the annihilation and creation
operators of the harmonic oscillator. The multi-mode Hamiltonian is
obtained summing over the number k of modes <cit.>, and
has the expression


    H=ω_AJ_z+∑_=1^kΩ_a_^†a_-1/√(N)∑_=1^kγ_(J_-+J_+)(a_+a_^†).



The k modes of the electromagnetic field are described in terms
of annihilation and creation operators for each mode a_,
a_^†, acting on the tensor product of k copies
of the Fock space 
and satisfying the commutation relations


    [a_,a_^†]=δ_,    [a_,a_]=[a_^†,a_^†]=0.



A two-level atom is described using the 1/2-spin matrices
S_z=1/2σ_z, S_±=1/2(σ_x± iσ_y)
(σ_x, σ_y and σ_z being the Pauli matrices),
which act on a two-dimensional complex Hilbert space ℂ^2
and satisfy the commutation relations


    [S_+,S_-]=2S_z,    [S_z,S_±]=± S_±.



When considering a system of N two-level atoms, we use the collective
spin operators J_z, J_-, J_+ defined as


    J_♢=S_♢⊗ I_2^⊗(N-1)+I_2⊗ S_♢⊗ I_2^⊗(N-2)
    +⋯+I_2^⊗(N-2)⊗ S_♢⊗ I_2+I_2^⊗(N-1)⊗ S_♢
where
I_2 is the identity operator on ℂ^2 and ♢∈{ z,-,+}.
These collective spin operators satisfy the commutation relations


    [J_+,J_-]=2J_z,    [J_z,J_±]=± J_±

and act, in principle, on the complex Hilbert space (ℂ^2)^⊗ N;
however, working with this space is physically equivalent to studying
a system of N fully distinguishable atoms, which we don't usually
have in the experimental setups used in the study of the QPT in the
Dicke model. To overcome this issue, we must use the common set of
eigenvectors {|j,m⟩} of the two
commuting observables J_z and J^2=1/2(J_+J_-+J_-J_+)+J_z^2,
where the label j is limited to the values j∈{ r,r+1,…,N/2}
(r=0 for even N and r=1/2 for odd N) and the label m∈ℤ
is constricted by |m|≤ j. These vectors do not form
a basis of (ℂ^2)^⊗ N for N>2,
as the dimension of their linear span is


    dim{span{{|j,m⟩} _|m|≤ j^j=r,…,N/2}} =∑_j=r^N/2(2j+1)≤2^N.









We will denote by ℋ_A the subspace of (ℂ^2)^⊗ N
generated by the states {|j,m⟩} _|m|≤ j^j=r,…,N/2.

There are two main results concerning the states {|j,m⟩} _|m|≤ j^j=r,…,N/2
and the space ℋ_A: the first comes from noticing that
[H,J^2]=0, which means that the label j of the
eigenvalues of J^2 remains constant during the system's evolution;
the second is the decomposition ℋ_A=,
where each ℋ_j is the subspace of dimension dim{ℋ_j} =2j+1
generated by the states {|j,m⟩} _|m|≤ j
with a fixed j. In this treatment, in order to study indistinguishable
atoms, we are ignoring the multiplicities g(j) of the irreducible
representations of SU(2), i.e. the number of times that
each ℋ_j appears in the full decomposition (ℂ^2)^⊗ N=.

To make it clear that the space ℋ_A is the one we must
work with when indistinguishable atoms are considered we should inquire
into the physical interpretation of the labels j and m. In order
to give a physical interpretation to the label j we must notice
that the energy of the atomic system is bounded by ± jω
independently of the number of atoms N (but with the restriction
j≤N/2), this leads us to interpret the quantity 2j
as the effective number of atoms in the system and define it as the
cooperation number. To make the notion of the cooperation
number more intuitive, Dicke, in his original paper <cit.>,
compares a state with j=0, which exists only for an even number
of atoms, with a classical system of an even number of oscillators
swinging in pairs oppositely phased. The interpretation of the label
m is clear from the definition of J_z: m=1/2(n_e-n_g),
where n_e and n_g are the number of atoms in the excited
and ground states, respectively.

In this paper we restrict our analysis to the space ℋ_A,
as it allows us to choose j as an initial condition (which will
remain constant) and work in ℋ_j, where the atoms are
indistinguishable.










§ METHODOLOGY


There have been various contributions to the study of the phase transition
in the Dicke model (and other two-level models) <cit.>
and different approaches such as Husimi function analysis <cit.>,
entropic uncertainty relations <cit.> and energy surface minimization
<cit.>, have been used
for its investigation.

In this work we use the energy surface minimization method, which
consists on minimizing the surface that is obtained by taking the
expectation value of the modeling Hamiltonian with respect to some
trial variational state. The strength of this method lies on the choice
of the trial state, as it is the latter, after minimization, the one
that will be modeling the ground state of the system.

Here we take a variational approach for both matter and radiation
fields, and show how to calculate the QPT of the system modeled by
the Hamiltonian H given in eq. (<ref>) via four means:


  * Using a tensor product of Heisenberg-Weyl HW(1) coherent states for
each mode of the electromagnetic field and SU(2) coherent states for
the atomic field as trial states, and analytically minimizing the
obtained energy surface with respect to its parameters.

  * Using a projection operator on HW(1) coherent states and SU(2) coherent
states to obtain trial states that preserve the parity symmetry of
the Hamiltonian with respect to the total excitation number of the
system (symmetry adapted states), and numerically minimize
the obtained energy surface with respect to its parameters.

  * Using symmetry adapted states, as in (2) above, to obtain the energy
surface and “minimize” it with the minimizing parameters obtained
in (1) above, thus allowing us to have analytic expressions for the
ground state.

  * Numerically diagonalizing the Hamiltonian, which gives us the exact
quantum solution.







  §.§.§ Coherent states (CS)


For each mode of the electromagnetic field the annihilation and creation
operators a_ and a_^†, appearing in
the modeling Hamiltonian H, satisfy the commutation relations (<ref>)
of the Lie algebra generators of the Heisenberg-Weyl group HW(1);
hence, a natural choice of a trial state for the radiation field is
a tensor product of k (number of modes) coherent states of HW(1)


    |α̅⟩ :=|α_1⟩⊗⋯⊗|α_k⟩ ,

where each |α_⟩ is defined as


    |α_⟩ :=e^α_a_^†-α_^*a_|0_⟩ =e^-|α_|^2/2∑_ν_=0^∞α_^ν_/√(ν_!)|ν_⟩ .



Furthermore, the commutation relations of the collective spin operators
J_-, J_+ and J_z (<ref>) are the same as the
ones of the Lie algebra generators of the special unitary group SU(2).
Thus, analogously as for the radiation field, we use the coherent
states of SU(2)


    |ξ⟩ _j:=|υtan|υ|/|υ|⟩ _j:=e^υ J_+-υ^*J_-|j,0⟩
    =1/(1+|ξ|^2)^j∑_m=0^2j2jm^1/2ξ^m|j,m-j⟩ .
as
trial states for the matter field.




  §.§.§ Symmetry adapted states (SAS)


The modeling Hamiltonian we are considering has a parity symmetry
given by [e^iπΛ,H]=0, where Λ=√(J^2+1/4)-1/2+J_z+∑_=1^ka_^†a_
is the excitation number operator with eigenvalues λ=j+m+∑_=1^kν_.
This symmetry allows us to classify the eigenstates of H in terms
of the parity of the eigenvalues λ; however, as states with
opposite symmetry are strongly mixed by the CS defined in the previous
section, we should then adapt this symmetry to the CS by projecting
them with the operator P_±=1/2(I± e^iπΛ),
i.e.


    |α̅,ξ_j⟩ _±:=𝒩_±P_±|α̅⟩⊗|ξ⟩ _j
    =𝒩_±(|α̅⟩⊗|ξ⟩ _j±|-α̅⟩⊗|-ξ⟩ _j),
with
𝒩_±=(2±2E(-cosθ)^2j)^-1/2
the normalization factors for the even (+) and odd (-) states (where
E=exp{ -2∑_=1^k|α_^2|}).

As we are interested in the ground state of the system, which has
an even parity, we only focus on the state |α̅,ξ_j⟩ _+.




  §.§.§ Entropy of entanglement (S_ε)


Entropy of entanglement is defined for a bipartite system as the Von
Neumann entropy of either of its reduced states, that is, if ρ
is the density matrix of a system in a Hilbert space ℋ=ℋ_1⊗ℋ_2,
its entropy of entanglement is defined as


    S_ε:=-Tr{ρ_1logρ_1} =-Tr{ρ_2logρ_2} ,

where ρ_1=Tr_2{ρ} and ρ_2=Tr_1{ρ}.

Our Hamiltonian H models a bipartite system formed by matter and
radiation subsystems, which means that their entropy of entanglement
can be used to see the influence of the QPT on its behavior; this
we do below.




  §.§.§ Fidelity between neighboring states (F)


Fidelity is a measure of the distance
between two quantum states; given |ϕ⟩ and
|φ⟩ it is defined as


    F(ϕ,φ):=|⟨ϕ|φ⟩|^2.



Across a QPT the ground state of a system suffers a sudden, drastic
change, thus it is natural to expect a drop in the fidelity between
neighboring states near the transition. This drop has been, in fact,
already shown to happen <cit.> for the case 2j=N.
We study it here also, and its behavior with the cooperation number.










§ RESULTS


Writing the complex labels α_ and ξ as α_=q_+ip_ι
and ξ=tan(θ/2)e^iϕ, with q_,p_∈ℝ,
θ∈[0,π), ϕ∈[0,2π), the
CS's energy surface is obtained by taking the expectation value of
the modeling Hamiltonian H with respect to the state |α̅⟩⊗|ξ⟩ _j,
and has the form


    ℋ_j,CS(q_,p_,θ,ϕ):=⟨α̅|⊗⟨ξ|_jH|α̅⟩⊗|ξ⟩ _j
    =-jω_Acosθ+∑_=1^kΩ_(q_^2+p_^2)-4j/√(N)sinθcosϕ∑_=1^kγ_q_.


The critical points which minimize it are then found to be


    θ_c=q__c=p__c=0,     ω_A≥8j/N∑_=1^kγ_^2/Ω_,




    .[ cosθ_c=Nω_A/8j(∑_=1^kγ_^2/Ω_)^-1,;                                  ;                          ϕ_c=0,π,;                                  ;     q__c=2jγ_/Ω_√(N)cosϕ_csinθ_c,;                                  ;                            p__c=0 ]}    ω_A<8j/N∑_=1^kγ_^2/Ω_.









Substituting these values into (<ref>) we obtain the energy
of the coherent ground state as a function of the Hamiltonian parameters,


    ℰ_CS(ω_A,γ_)={ -jω_A    ,     δ≥1
     -jω_A/2(1/δ+δ)    ,     δ<1,
    .

where we have defined δ=Nω_A/8jς
with ς=∑_=1^kγ_^2/Ω_.
Using the information of this coherent ground state we also obtain
the expectation values of the atomic relative population operator
J_z and of the number of photons of mode  operator ν_:=a_^†a_:


    ⟨ J_z⟩ _CS(ω_A,γ_)={ -j    ,     δ≥1
     -jδ    ,     δ<1,
    .




    ⟨ν_⟩ _CS(ω_A,γ_,γ_)={ 0    ,     δ≥1
    γ_^2/Ω_^2jω_A/2ς(1/δ-δ)    ,     δ<1.
    .



Analogously as for the CS's energy surface, the SAS's energy surface
is obtained by taking the expectation value of the modeling Hamiltonian
H with respect to the state |α̅,ξ_j⟩ _+,
and has the more complicated form


    ℋ_j,SAS(q_,p_,θ,ϕ):=⟨α̅,ξ_j|_+H|α̅,ξ_j⟩ _+
    =(1+E(-cosθ)^2j-2/1+E(-cosθ)^2j)(-jω_Acosθ)
    +(1-E(-cosθ)^2j/1+E(-cosθ)^2j)∑_ℓ=1^kΩ_ℓ(q_ℓ^2+p_ℓ^2)
    -4j/√(N)sinθ∑_ℓ=1^k{cosϕγ_ℓq_ℓ+E(-cosθ)^2j-1sinϕγ_ℓp_ℓ/1+E(-cosθ)^2j}.


As a first approximation, we may substitute the critical values obtained
for the CS's energy surface into (<ref>), we obtain the trial
state which approximates the lowest symmetry-adapted energy state,
and with respect to which we evaluate the expectation values of the
observables H, J_z and ν_:


    ℰ_SAS(ω_A,γ_)
    ={[                             -jω_A                            ,  δ≥1; ; -jω_A[δ(1+ε(-δ)^2j-2/1+ε(-δ)^2j).;                     .+1/2(1/δ-δ)]                           ,  δ<1, ].



    ⟨ J_z⟩ _SAS(ω_A,γ_)
    ={ -j    ,     δ≥1
     -jδ(1+ε(-δ)^2j-2/1+ε(-δ)^2j)    ,     δ<1,.



    ⟨ν_⟩ _SAS(ω_A,γ_,γ_)
    ={ 0    ,     δ≥1
    γ_^2/Ω_^2jω_A/2ς(1/δ-δ)(1-ε(-δ)^2j/1+ε(-δ)^2j)    ,     δ<1,.
where
ε=exp{-jω_Aσ/ς}
with σ=∑_=1^kγ_^2/Ω_^2.

Of course, we can minimize eq. (<ref>) numerically for the
SAS and obtain the expectation value of the relevant matter and field
observables.

In our numerical analysis we study the case with two modes of the
radiation field, as it is the maximum number of orthogonal modes that
can be present in a 3D cavity with the restrictions that the modes
interact with the electric dipole moment of the atoms and to be in
resonance with the frequency associated with the energy difference
between the two levels of the atoms. This latter restriction is just
considered to have the maximum transition probability between states.

For the exact quantum solution we must resort to numerical diagonalization
of the Hamiltonian and use the lowest eigenstate to compute the expectation
values of the relevant observables.

The results, properties of the ground state related to the CS, those
related to the SAS using the critical points of the CS (which have
the advantage of also providing analytical solutions), those of the
SAS minimized numerically and those of the quantum solution through
numerical diagonalization, are shown in figures <ref> - <ref>
and are discussed below.

One advantage of having analytical solutions is, of course, that the
order of the transition may be easily found. Equations (<ref>)
and (<ref>) show a second-order QPT at δ=Nω_A/8jς=1
with the CS and SAS using CS's minima (SASc) approximations. In figure
<ref> it can be seen that the data of the SAS using numerical
minimization (SASn) has a small discontinuity (the QPT) at γ_2≈1.485
for j=5 and γ_2≈1.015 for j=9, while in figure
<ref> this discontinuity is at ω_A≈0.975 for
j=5 and ω_A≈1.965 for j=9. Note that the SASn
solution always approximate better the exact quantum result, as the
cooperation number increases this approximation gets better, in fact,
for 2j=18=N the loci of the separatrix between the normal and collective
regions for the quantum and SASn solutions are indistinguishable (except
in the zoomed inset). The true loci of the QPT may be found through
the fidelity: figures <ref> and <ref> show the fidelity
between neighboring states of the quantum solution, where the exact
QPT is characterized by the minimum, which is localized at γ_2≈1.550
for j=5, γ_2≈1.031 for j=9 in figure <ref>;
and ω_A≈0.817 for j=5, ω_A≈1.870
for j=9 in figure <ref>.

The discrepancies between the transition values of the SASc approximation
and the exact quantum solution become obvious when looking at figures
<ref> and <ref>, where the fidelity between SASc and
the quantum solution drops (and oscillates) in a vicinity of the separatrix.
Therefore, we conclude that SASc offer a good approximation (with
an analytic expression) to the exact quantum solution far from the
QPT for low cooperation numbers, but as j→∞, the
interval where the SASc fail to reproduce the correct behavior, becomes
smaller.

Figures <ref> and <ref> show the fidelity drop at
the separatrix for the SASn. The resemblance to figures <ref>
and <ref> is uncanny, showing the benefits of restoring the
Hamiltonian symmetry into the trial variational states. This improvement
comes with the disadvantage of losing the analytic expression, but
still has an advantage over the quantum solution: the computational
time. SASn are obtained by numerically minimizing a real function,
which is far easier to do (computationally speaking) than numerically
diagonalizing the Hamiltonian matrix.







Figures <ref> - <ref> show the comparison between the
different approximations to the ground state: CS, SASc, SASn and quantum
solution. We show the behavior of ℰ:=⟨ H⟩,
⟨ J_z⟩ and ⟨ν_⟩
as functions of the atomic frequency ω_A and one of the
coupling constants γ_2, for different cooperation numbers.
It can be noticed that the discontinuity in the second derivative
of the energy (as modeled with CS and SASc) translates into a discontinuity
in the first derivative of ⟨ J_z⟩ and
⟨ν_⟩, thus characterizing the
QPT by means of an abrupt change in the expectation values of the
observables. In general, it can be observed that the four methods
(CS, SASc, SASn and quantum solution) converge in the limit δ→0,
where the case j→∞ is particularly interesting as
the interval around the QPT, where all the approximations fail to
reproduce the correct behavior, becomes smaller.

It is worth mentioning the significance and importance of figures
<ref> and <ref> as they show aspects of the multi-mode
Dicke model which are not present in the single-mode case. In figure
<ref> it is shown how the different modes of radiation (orthogonal
in principle) interact through the matter field, analogously as it
occurs with different atoms interacting through the radiation field.
On the other hand, figure <ref> shows (pictorically) the phase
diagrams of the two-mode system, in which it can be observed that
any two points in the super-radiant region can be joined by a trajectory
that does not cross the normal region, a characteristic that the single-mode
system does not have. 

Figures <ref> and <ref> show the comparison between
SASc, SASn and the quantum solution for the entropy of entanglement
S_ε as a function of the atomic frequency ω_A
and one of the coupling constants γ_2, using different cooperation
numbers. A characterization of the QPT can be made by observing that
the entropy of entanglement obtained using the quantum solution shows
a maximum at the transition, an attribute that SASc and SASn approximations
fail to reproduce.















§ DISCUSSION AND CONCLUSIONS


From figures presented we conclude that SASc offer a good approximation
(with an analytic expression) to the exact quantum solution far from
the QPT for low cooperation numbers, but as j→∞,
the interval where the SASc fail to reproduce the correct behavior,
becomes smaller.

On the other hand, the SASn provide a better approximation to the
quantum solution. This improvement comes with the disadvantage of
losing the analytic expression, but still has the advantage over the
quantum solution of the computational time. SASn are obtained by numerically
minimizing a real function, which is far easier to do (computationally
speaking) than numerically diagonalizing the Hamiltonian matrix.

A characterization of the QPT can be made by looking at the entropy
of entanglement; that obtained using the quantum solution shows a
maximum at the transition, an attribute that SASc and SASn approximations
fail to reproduce.

The behavior of the expectation values of the relevant observables
of the system ⟨ H⟩, ⟨ J_z⟩
and ⟨ν_⟩, is also affected by
the QPT (figures <ref> - <ref>), thus allowing us to
characterize the QPT by means of its influence over the observables.
In general, it can be observed that the four methods (CS, SASc, SASn
and quantum solution) converge in the limit δ→0,
where the case j→∞ is particularly interesting as
the interval around the QPT, where all the approximations are weaker,
becomes smaller.

In conclusion, we have shown how the use of variational states to
approximate the ground state of a system can be useful to characterize
the QPT in a multi-mode Dicke model using the expectation value of
the observables relevant to the system and the entropy of entanglement
between matter and radiation. We have also introduced a not very commonly
used dependence: the cooperation number, showing its influence over
the behavior of the system, paying particular attention to the QPT
and the accuracy of the used approximations. Some aspects of the multi-mode
Dicke model which are not present in the single-mode case were also
briefly discussed.



We thank R. López-Peña, O. Castaños and S. Cordero for their comments
and discussion. This work was partially supported by DGAPA-UNAM under
project IN101217. L. F. Q. thanks CONACyT-México for financial support
(Grant #379975).



10
	
key-1R. H. Dicke, http://journals.aps.org/pr/abstract/10.1103/PhysRev.93.99Phys. Rev. 93, 99 (1954).

key-2K. Hepp and E. H. Lieb, http://www.sciencedirect.com/science/article/pii/0003491673900390Ann. Phys. 76, 360 (1973).

key-3K. Hepp and E. H. Lieb, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.8.2517Phys. Rev. A 8, 2517 (1973).

key-4Y. K. Wang and F. T. Hioe, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.7.831Phys. Rev. A 7, 831 (1973).

key-5K. Baumann, C. Guerlin, F. Brennecke and T. Esslinger,
http://www.nature.com/nature/journal/v464/n7293/abs/nature09009.htmlNature (London) 464, 1301 (2010).

key-6D. Nagy, G. Kónya, G. Szirmai and P. Domokos, http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.104.130401Phys. Rev. Lett. 104, 130401 (2010).

key-7T. Brandes, http://www.sciencedirect.com/science/article/pii/S0370157304005496Phys. Rep. 408 , 315 (2005).

key-8G. Chen, Z. Chen and J. Liang, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.76.055803Phys. Rev. A 76 , 055803 (2007).

key-9N. Lambert, C. Emary and T. Brandes, http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.92.073602Phys. Rev. Lett. 92, 073602 (2004).

key-10N. Lambert, C. Emary and T. Brandes, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.71.053804Phys. Rev. A 71, 053804 (2005).

key-11E. Nahmad-Achar, S. Cordero, O. Castaños and R. López-Peña,
http://iopscience.iop.org/article/10.1088/0031-8949/90/7/074026Phys. Scr. 90, 074026 (2015).

key-12E. Nahmad-Achar, S. Cordero, R. López-Peña and O.
Castaños, http://iopscience.iop.org/article/10.1088/1751-8113/47/45/455301/meta;jsessionid=100FF83C0FD0D059102F9D188FAFA77E.ip-10-40-1-98J. Phys. A: Math. Theor. 47, 455301 (2014).

key-13S. Cordero, O. Castaños, R. López-Peña and E. Nahmad-Achar,
http://iopscience.iop.org/article/10.1088/1751-8113/46/50/505302J. Phys. A: Math. Theor. 46, 505302 (2013).

key-14J. Reslen, L. Quiroga and N. F. Johnson, http://iopscience.iop.org/article/10.1209/epl/i2004-10313-4/metaEurophys. Lett. 69, 1 (2005).

key-15P. Zanardi and N. Paunković, http://journals.aps.org/pre/abstract/10.1103/PhysRevE.74.031123Phys. Rev. E 74, 031123 (2006).

key-16H. Goto and K. Ichimura, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.77.053811Phys. Rev. A 77, 053811 (2008).

key-17C. Emary and T. Brandes, http://journals.aps.org/pre/abstract/10.1103/PhysRevE.67.066203Phys. Rev. E 67, 066203 (2003).

key-18C. Emary and T. Brandes, http://journals.aps.org/pre/abstract/10.1103/PhysRevE.67.066203Phys. Rev. Lett. 90, 044101 (2003).

key-19E. Romera, R. del Real and M. Calixto, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.85.053831Phys. Rev. A 85, 053831 (2012).

key-20E. Romera, M. Calixto and Á. Nagy, http://iopscience.iop.org/article/10.1209/0295-5075/97/20011/metaEurophys. Lett. 97, 20011 (2012).

key-21O. Castaños, R. López-Peña, E. Nahmad-Achar, J. G.
Hirsch, E. López-Moreno and J. E. Vitela, http://iopscience.iop.org/article/10.1088/0031-8949/79/06/065405/metaPhys. Scr. 79, 065405 (2009).

key-22O. Castaños, E. Nahmad-Achar, R. López-Peña and J.
G. Hirsch, http://iopscience.iop.org/article/10.1088/0031-8949/80/05/055401/metaPhys. Scr. 80, 055401 (2009).

key-23O. Castaños, E. Nahmad-Achar, R. López-Peña, and
J. G. Hirsch, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.83.051601Phys. Rev. A 83, 051601 (2011).

key-24O. Castaños, E. Nahmad-Achar, R. López-Peña, and
J. G. Hirsch, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.84.013819Phys. Rev. A 84, 013819 (2011).

key-25O. Castaños, E. Nahmad-Achar, R. López-Peña, and
J. G. Hirsch, http://journals.aps.org/pra/abstract/10.1103/PhysRevA.86.023814Phys. Rev. A 86, 023814 (2012).

key-26J. G. Hirsch, O. Castaños, E. Nahmad-Achar and R.
López-Peña, http://iopscience.iop.org/article/10.1088/0031-8949/2013/T153/014033/metaPhys. Scr. 2013, 014033 (2013).

