text:
# 1 Introduction

## 1 Introduction

Resonant interaction of an alternating electromagnetic field with an atom makes it possible to select a pair of atomic levels and to treat the atom as an effectively two-level system . This resonant interaction is in the basis of one of the main directions of optics. The possibility of realizing atomic resonance is due to the existence of discrete energy levels of electrons in an atom, the energy spectrum being described by the Schrödinger equation. The Bose-condensed atoms are also described by the Schrödinger equation, though nonlinear, which is often termed the Gross-Pitaevskii equation \[2–4\]. Trapped atoms, because of a confining potential, also possess a discrete spectrum . A specific feature is that this is a spectrum of collective nonlinear states of an ensemble of coherent atoms. The collective coherent and nonlinear nature of these states is what makes them principally different from single-particle linear states of an individual atom. To emphasize this difference, the collective states of Bose-condensed atomic gases are called nonlinear coherent modes . Since the spectrum of these modes, pertaining to trapped atoms, is discrete, it has been suggested  that a resonant coupling of two such modes is feasible by means of an alternating field with a frequency tuned to the transition frequency between the chosen levels. This type of resonance for exciting nonlinear coherent modes can be named coherent resonance. This phenomenon, if realised, could open a whole new field of possible applications, analogous to those accomplished with atomic resonance.

In the present communication, we investigate the conditions required for realizing the coherent resonance. We propose a new method for stabilizing Bose-condensed gas with attractive interactions by transferring atoms to an excited coherent mode. We study dynamic critical phenomena, discovered numerically \[6–8\], and give their explanation. We also describe some new properties of Bose condensate subject to coherent resonance, such as interference effects, atomic squeezing, and multiparticle coherent entanglement.

## 2 Coherent Resonance

We consider a system of trapped Bose-condensed atoms at low temperature, when the Bose gas in a coherent state. The coherent-field wave function satisfies the Gross-Pitaevskii equation

$$i\mathrm{}\frac{}{t}\phi (𝐫,t)=\left(\widehat{H}[\phi ]+\widehat{V}\right)\phi (𝐫,t),$$
(1)
in which the coherent field is normalized to unity, $`\phi ^2=1`$. The nonlinear Hamiltonian is

$$\widehat{H}[\phi ]\frac{\mathrm{}^2^2}{2m_0}+U(𝐫)+NA_s|\phi (𝐫,t)|^2,$$
(2)
where $`U(𝐫)`$ is a trapping potential; $`m_0`$ is atomic mass; $`N`$ is the number of atoms; $`A_s4\pi \mathrm{}^2a_s/m_0`$, with $`a_s`$ being a scattering length. The alternating field is

$$V(𝐫,t)=V_1(𝐫)\mathrm{cos}\omega t+V_2(𝐫)\mathrm{sin}\omega t.$$
(3)
Note that Eq. (1) is an exact equation for the coherent field , and $`\phi (𝐫,t)`$ does not need to be interpreted as an average of a field operator, for which case Eq. (1) would be only a mean-field approximation.

The nonlinear coherent modes are the stationary solutions to the eigenproblem

$$\widehat{H}[\phi _n]\phi _n(𝐫)=E_n\phi _n(𝐫).$$
(4)
Since the Hamiltonian (2) is nonlinear, the modes are not necessary orthogonal to each other, so that the scalar product

$$(\phi _m,\phi _n)\phi _m^{}(𝐫)\phi _n(𝐫)𝑑𝐫$$

is not compulsory a Kroneker delta. But each mode can always be normalized, $`\phi _n^2=1`$. Some properties of the nonlinear modes have been discussed \[5–8,10–13\] and dipole mode has been observed experimentally .

Let us select two spectrum levels, with the energies $`E_1`$ and $`E_2`$, such that $`E_1<E_2`$. The related transition frequency is

$$\omega _{21}\frac{1}{\mathrm{}}(E_2E_1).$$
(5)
The frequency of the alternating field (3) is tuned close to this transition frequency (5), which implies the resonance condition

$$\left|\frac{\mathrm{\Delta }\omega }{\omega }\right|1,\mathrm{\Delta }\omega \omega \omega _{21}.$$
(6)
The solution to the evolution equation (1) can be presented in the form

$$\phi (𝐫,t)=\underset{n}{}c_n(t)\phi _n(𝐫)\mathrm{exp}\left(\frac{i}{\mathrm{}}E_nt\right).$$
(7)
We assume that $`c_n(t)`$ is slow as compared to the exponential, such that

$$\frac{\mathrm{}}{E_n}\left|\frac{dc_n}{dt}\right|1.$$
(8)
There are two transition amplitudes in the system, one is the internal transition amplitude

$$\alpha _{mn}A_s\frac{N}{\mathrm{}}|\phi _m(𝐫)|^2\left[2|\phi _n(𝐫)|^2|\phi _m(𝐫)|^2\right]𝑑𝐫,$$
(9)
caused by atomic interactions, and another is the external transition amplitude

$$\beta _{mn}\frac{1}{\mathrm{}}\phi _m^{}(𝐫)\left[V_1(𝐫)iV_2(𝐫)\right]\phi _n(𝐫)𝑑𝐫,$$
(10)
due to the driving resonant field. Similarly to atomic resonance, in order to avoid power broadening, the transition aplitudes are to be smaller than the transition frequencies,

$$\left|\frac{\alpha _{mn}}{\omega _{mn}}\right|1,\left|\frac{\beta _{mn}}{\omega _{mn}}\right|1.$$
(11)
Under conditions (6), (8), and (11), the Gross-Pitaevskii equation (1) reduces to the set of equations

$$i\frac{dc_1}{dt}=\alpha _{12}|c_2|^2c_1+\frac{1}{2}\beta _{12}c_2e^{i\mathrm{\Delta }\omega t},i\frac{dc_2}{dt}=\alpha _{21}|c_1|^2c_2+\frac{1}{2}\beta _{12}^{}c_1e^{i\mathrm{\Delta }\omega t}$$
(12)
for an effective two-mode system, where $`|c_1|^2+|c_2|^2=1`$. This reduction is valid provided all transition amplitudes, involved in Eq. (12), are small in the sense of condition (11), which yields the inequalities

$$\left|\frac{\alpha _{12}}{\omega _{21}}\right|1,\left|\frac{\alpha _{21}}{\omega _{21}}\right|1,\left|\frac{\beta _{12}}{\omega _{21}}\right|1.$$
(13)
To analyse the above inequalities, we consider a cylindrical trap modelled by the harmonic trapping potential

$$U(𝐫)=\frac{m_0}{2}\left(\omega _r^2r_x^2+\omega _r^2r_y^2+\omega _z^2r_z^2\right),$$

with the aspect ratio

$$\nu \omega _z/\omega _r.$$
(14)
Introduce the dimensionless coupling parameter

$$g4\pi \frac{a_s}{l_r}N\left(l_r\sqrt{\frac{\mathrm{}}{m_0\omega _r}}\right).$$
(15)
The spectrum of nonlinear coherent modes is defined by the eigenproblem (4). The external transition amplitude (10) can always be made sufficiently small by regulating the amplitude of the resonant field. We need to calculate the internal transition amplitude (9) and to compare it with a chosen transition frequency (5). Calculations can be accomplished by means of the optimized perturbation theory \[15–17\] (for reviews see ). The first two inequalities (13) can be valid only outside the radius of convergence of the strong-coupling expansion. This imposes the restriction

$$|g\nu |<\frac{[2p^2+(q\nu )^2]^{5/4}}{14p\sqrt{q}I_{nmj}}$$
(16)
on the coupling parameter (15). Here the integral $`I_{nmj}(|\psi _{nmj}|^2,|\psi _{nmj}|^2)`$, and $`p2n+|m|+1`$, $`q2j+1`$ are the combinations of the radial, $`n`$, azimuthal, $`m`$, and axial, $`j`$, quantum numbers. The restriction (16) can be rewritten as the limitation on the admissible number of particles $`N<N_0`$, with the limiting number

$$N_0=\frac{[2p^2+(q\nu )^2]^{5/4}}{56\pi \nu p\sqrt{q}I_{nmj}}\left|\frac{l_r}{a_s}\right|.$$
(17)
For the ground-state level, this reduces to

$$N_0=\sqrt{\frac{\pi }{2}}\frac{(2+\nu ^2)^{5/4}}{14\nu }\left|\frac{l_r}{a_s}\right|.$$
(18)
Thus, the number of particles, which can be resonantly transferred to an excited coherent mode, is limited by a number $`N_0`$, which depends on the quantum numbers of the coupled modes. The number $`N_0`$ also essentially depends on the trap shape through the aspect ratio (14).

It turns out that the limiting numbers (17) and (18) are close to the critical values defining the maximal number of atoms with a negative scattering length, which can be condensed in a trap. Examples of atoms with attractive interactions are 7Li (see review ) and 85Rb (see Ref. ). Such atoms can form condensates only if their number does not exceed a critical value \[5,22–25\]. In that case, the condensate forms a long-lived metastable state, with a very slow decay, caused by quantum tunneling \[26–28\], with the tunneling probability being negligible as compared to the probability of escaping from the trap because of depolarizing collisions . But if the number of atoms is larger than critical, the condensate collapses. With a supply of atoms from an external source, say by a continuous loading from another trap , the condensate grows again and, after surpassing the critical number, again collapses. Thus a series of growths and collapses take place, as was observed experimentally  in Bose-condensed 7Li and described theoretically by means of the Gross-Pitaevskii equation supplemented by relaxation terms and by the related rate equations. Explosion of an attractive condensate of 85Rb atoms by manipulating the atomic interactions with an external magnetic field near a Feshbach resonance  has also been observed .

To increase the number of atoms in a condensate with attractive interactions, it was suggested  to drive a quadrupole collective excitation. We may notice that transferring atoms to an excited nonlinear coherent mode may also stabilize such a condensate. This follows from the form (17), which shows that the limiting number $`N_0`$ increases for higher excited modes as

$$N_0(2n+|m|+1)^{3/2},N_0(2j+1)^2.$$

To estimate the number of atoms for highly excited modes, one could also use an optimized quasiclassical approximation . As Eq. (18) demonstrates, the limiting number $`N_0`$ is larger for essentially anisotropic traps, that is, for cigar-shape $`(\nu 1)`$ and disk-shape $`(\nu 1)`$ traps.

The temporal behaviour of the system under coherent resonance is described by the evolution equations (12). Numerical solution of these equations, for the varying amplitude of the resonant field, displays dynamical critical phenomena \[6–8\]. Here we present a general picture of these phenomena and elucidate their origin.

Let us present the population amplitudes in the form

$$c_1=\sqrt{\frac{1s}{2}}\mathrm{exp}\left\{i\left(\pi _1+\frac{\mathrm{\Delta }\omega }{2}t\right)\right\},c_2=\sqrt{\frac{1+s}{2}}\mathrm{exp}\left\{i\left(\pi _2\frac{\mathrm{\Delta }\omega }{2}t\right)\right\},$$
(19)
where $`\pi _1=\pi _1(t)`$ and $`\pi _2=\pi _2(t)`$ are the phases. The variable

$$s|c_2|^2|c_1|^2$$
(20)
is the population difference. Introduce the notation for the average internal transition amplitude

$$\alpha \frac{1}{2}(\alpha _{12}+\alpha _{21})$$
(21)
and let us take into account that the external transition amplitude (10) is, in general, complex valued,

$$\beta _{12}=\beta e^{i\gamma },\beta |\beta _{12}|.$$
(22)
Also, define the effective detuning

$$\delta \mathrm{\Delta }\omega +\frac{1}{2}(\alpha _{12}\alpha _{21}).$$
(23)
And introduce the phase difference

$$x\pi _2\pi _1+\gamma .$$
(24)
Substituting Eqs. (19) into the evolution equations (12), we come to the set of equations for the population difference (20) and phase difference (24),

$$\frac{ds}{dt}=\beta \sqrt{1s^2}\mathrm{sin}x,\frac{dx}{dt}=\alpha s+\frac{\beta s}{\sqrt{1s^2}}\mathrm{cos}x+\delta .$$
(25)
Numerical solution of Eqs. (25) demonstrates the existence of dynamical critical phenomena, when a tiny variation of parameters results in sharp changes of temporal behaviour of $`s(t)`$ and $`x(t)`$. We illustrate this in the following figures, where time is measured in units of $`\alpha ^1`$. There are two independent dimensionless parameters

$$b\frac{\beta }{\alpha },\epsilon \frac{\delta }{\alpha },$$
(26)
which can be varied. In what follows, we vary $`b`$ in the interval $`[1,1]`$ and keep $`|\epsilon |1`$. And there are also two initial conditions $`s_0=s(0)`$ and $`x_0=x(0)`$.

In Fig. 1, we set $`s_0=1`$, $`x_0=0`$, and $`\epsilon =0`$, while varying the parameter $`b`$ corresponding to the amplitude of the pumping resonant field. In the range $`0<b<0.5`$, the population difference oscillates between $`1`$ and zero and the phase difference monotonically decrease, as is shown in Fig. 1a and 1b. The critical point $`b_c=0.5`$ separates the regions of two qualitatively different temporal behaviours. After $`b`$ surpasses $`b_c`$, the population difference starts oscillating between $`1`$ and $`+1`$, while the phase difference becomes a periodic function, as is demonstrated in Figs. 1c to 1e.

As follows from the definition of the internal transition amplitude (9), its value may become negative, either for strongly energetically separated modes or for atoms with attractive interactions, when $`A_s<0`$. In such a case, the parameter (21) is negative, $`\alpha <0`$, and so is the parameter $`b<0`$ from Eq. (26). The variation of $`b`$ in the negative region is also accompanied by dynamic critical phenomena, as for positive $`b`$, though the behaviour of the phase difference is not the same as for $`b>0`$. This is illustrated by Fig. 2, where we set the same initial conditions $`s_0=1`$ and $`x_0=0`$, and zero detuning. For $`b`$ in the interval $`0.5<b<0`$, the population difference is in the range $`1s0`$, as is shown in Figs. 2a and 2b. The critical point now is $`b_c=0.5`$. After crossing the critical point, when $`b<b_c`$, the amplitude of the population difference increases by a jump, so that $`s(t)`$ oscillates now in the diapason $`1s1`$. The phase difference is an oscillating function for both $`b>b_c`$ and $`b<b_c`$, but its shape changes drastically when crossing $`b_c`$, as can be seen in Figs. 2c and 2d.

In this way, there are two critical lines for each given pair of initial conditions $`s_0`$ and $`x_0`$. Changing initial conditions changes the location of the critical line on the manifold of the parameters $`b`$ and $`\epsilon `$. Thus, in Fig. 3, we show the course of events for the initial conditions $`s_0=0.8`$ and $`x_0=0`$, with the detuning $`\epsilon =0`$ and for the parameter $`b`$ varying in the whole range $`1b1`$. Under the fixed detuning $`\epsilon =0`$ and given initial conditions, there are two critical points $`b_{c1}=0.8`$ and $`b_{c2}=0.2`$.

For a nonzero detuning, the location of the critical line changes, but the overall picture remains the same. A finite detuning also makes the shape of the function $`x(t)`$ a little asymmetric, as is shown in Fig. 4.

To understand the origin of these critical phenomena, we studied the phase portrait on the plane $`\{s,x\}`$ for different parameters $`b`$ and $`\epsilon `$. Examples are demonstrated in Fig. 5. We also accomplished the stability analysis of fixed points. This investigation clarify the origin of the dynamic critical phenomena. The latter are related to the existence of the saddle separatrix defined by the equation

$$\frac{s^2}{2}b\sqrt{1s^2}\mathrm{cos}x+\epsilon s=b.$$
(27)
The separatrices divide the phase plane $`\{s,x\}`$ into the regions with different dynamics. As is evident from Eq. (27), changing any of the parameters $`b`$ or $`\epsilon `$ shifts the separatrices. When a separatrix crosses an initial point $`\{s_0,x_0\}`$, the system trajectory passes to another region of the phase plane, as a result of which the system dynamics changes qualitatively. The separatrix crossing of an initial point corresponds to the critical line

$$\frac{s_0^2}{2}b\sqrt{1s_0^2}\mathrm{cos}x_0+\epsilon s_0=|b|$$
(28)
on the manifold of the parameters $`b`$ and $`\epsilon `$. It is this separatrix crossing effect that causes the appearance of the dynamic critical phenomena. It can be shown  that a time-averaged system displays critical phenomena, typical of phase transitions in stationary statistical systems, on the same critical line.

There exist several more interesting effects occurring in the process of coherent resonance. Not to overload the present communication, we mention these effects only briefly. A more detailed consideration will be done in separate publications.

### 2.1 Interference effects

Coherent resonance couples two selected nonlinear modes, leaving other modes unpopulated. Therefore, the coherent field (7) is, effectively, a sum

$$\phi (𝐫,t)=\phi _1(𝐫,t)+\phi _2(𝐫,t)$$
(29)
of two terms

$$\phi _j(𝐫,t)=c_j(t)\phi _j(𝐫)\mathrm{exp}\left(\frac{i}{\mathrm{}}E_jt\right)(j=1,2).$$

Using this, one can define the interference pattern

$$\rho _{int}(𝐫,t)\rho (𝐫,t)\rho _1(𝐫,t)\rho _2(𝐫,t),$$
(30)
in which

$$\rho (𝐫,t)=|\phi (𝐫,t)|^2,\rho _i(𝐫,t)=|\phi _i(𝐫,t)|^2,$$

and the interference current

$$𝐣_{int}(𝐫,t)𝐣(𝐫,t)𝐣_1(𝐫,t)𝐣_2(𝐫,t),$$
(31)
where

$$𝐣(𝐫,t)=\frac{\mathrm{}}{m_0}\mathrm{Im}\phi ^{}(𝐫,t)\phi (𝐫,t),𝐣_i(𝐫,t)=\frac{\mathrm{}}{m_0}\mathrm{Im}\phi _i^{}(𝐫,t)\phi _i(𝐫,t).$$

These quantities show that the total density of atoms in a trap is not simply a sum of the partial mode densities but it includes also interference fringes and that, because of the different mode topology, there arises a local interference current. The interference effects can be observed experimentally.

### 2.2 Atomic squeezing

The evolution equations for the effective two-mode system, obtained by means of the coherent resonance, can be interpreted as the equations for the statistical averages of quasispin operators

$$S_\alpha \underset{i=1}{\overset{N}{}}S_i^\alpha ,S_\pm S_x\pm iS_y,$$

where $`\alpha =x,y,z`$ and $`S_i^\alpha `$ is an $`\alpha `$-component of the spin-1/2 operator. For instance, the population difference (20) can be written as the average

$$s=\frac{2}{N}<S_z>.$$

The squeezing factor may be defined  as

$$Q_z\frac{2\mathrm{\Delta }^2(S_z)}{\sqrt{<S_x>^2+<S_y>^2}}=\frac{2\mathrm{\Delta }^2(S_z)}{|<S_\pm >|},$$
(32)
where $`\mathrm{\Delta }^2(S_z)<S_z^2><S_z>^2`$. One says that $`S_z`$ is squeezed with respect to $`S_\pm `$ if $`Q_z<1`$. Calculating the factor (32), we find

$$Q_z=\sqrt{1s^2}.$$
(33)
Since $`s^21`$, one has $`Q_z1`$. The maximal squeezing is achieved when one of the modes is completely populated, i.e. $`s=\pm 1`$. Squeezed $`S_z`$ means that the population difference can be measured with a better accuracy than the phase difference.

### 2.3 Multiparticle entanglement

The system of interacting atoms possessing internal states may form multiparticle entangled states . In our case, atoms have not internal but collective coherent modes. This makes the principal difference of the resonant Bose condensate having collective nonlinear modes from atoms with internal single-particle states. But the resonant Bose condensate also forms a multiparticle entangled state. This follows from the Schrödinger representation for the wave function of $`N`$ atoms in the resonant condensate,

$$\mathrm{\Psi }_N(𝐫_1,𝐫_2,\mathrm{},𝐫_N,t)=c_1(t)\underset{i=1}{\overset{N}{}}\phi _1(𝐫_i)+c_2(t)\underset{i=1}{\overset{N}{}}\phi _2(𝐫_i),$$
(34)
where $`\phi _j(𝐫)`$ is a coherent $`j`$-mode, with $`j=1,2`$. The function (34) cannot be factorized into a product of single-particle functions, provided that $`c_10`$ and $`c_20`$. The maximal entanglement is reached for $`|c_1|=|c_2|=1/\sqrt{2}`$, that is, when $`s=0`$. Since the coefficients $`c_1(t)`$ and $`c_2(t)`$ are defined by the evolution equations (12) or (25), from where

$$|c_1|^2=\frac{1s}{2},|c_2|^2=\frac{1+s}{2},$$
(35)
their values can be manipulated by applying the resonant field and switching it off at the appropriate moment. In this way, it is possible to create any degree of entanglement and then one can disentangle the created entangled state. The resonant field, thus, acts as a quantum disentanglement eraser .

## 3 Conclusion

In trapped Bose condensates, the effect of coherent resonance is feasible. This is achieved by applying an alternating field with the frequency tuned to the transition frequency between two nonlinear modes. To preserve the resonance condition requires that the number of atoms be limited by a maximal value. Estimates show that this value is around $`10^3`$ for a typical spherical trap and for alkali atoms. The limiting number is higher for cigar-shape and disk-shape traps, and can be about $`10^5`$ atoms. This number also increases for higher excited nonlinear modes. Temporal behaviour of fractional populations exhibits dynamic critical phenomena occurring on a critical line in the parametric manifold. The origin of these critical phenomena is the saddle separatrix crossing by the starting point of a trajectory. Actually, there exists a whole bunch of critical lines, each of which corresponds to a different starting point.

As a result of coherent resonance in a system of Bose-condensed atoms, a new type of matter is formed, which can be called resonant Bose condensate. This is analogous to a resonant atom, with the principal difference that the resonant condensate is a coherent multiparticle system. Collective nature of the latter causes the appearance of interference patterns, interference current, atomic squeezing, and of multiparticle coherent entanglement.

We have considered here a single-component Bose gas. Possible extensions of the theory could be to binary mixtures of Bose condensates and, generally, to multicomponent condensates composed of different atomic species. Another possible generalization could be to the systems of cold atoms in optical lattices  and to radiating atoms interacting with electromagnetic field . We think that the effect of coherent resonance can find numerous applications, such as creating selected modes of atoms lasers, information processing, and quantum computing.

Acknowledgement

The work was accomplished in the Research Center for Optics and Photonics, University of São Paulo, São Carlos. Financial support from the São Paulo State Research Foundation is appreciated.

Figure Captions

Fig. 1. The population difference $`s(t)`$ (dashed line) and phase difference $`x(t)`$ (solid line) as functions of time, measured in units of $`\alpha ^1`$, for the zero detuning $`\epsilon =0`$, initial conditions $`s_0=1`$, $`x_0=0`$, and varying amplitude of the resonant field: (a) $`b=0.470`$; (b) $`b=0.495`$; (c) $`b=0.501`$; (d) $`b=0.7`$; (e) $`b=1`$.

Fig. 2. Temporal behaviour of the population difference $`s(t)`$ (dashed line) and phase difference $`x(t)`$ (solid line), under the same consitions as for Fig. 1, but for the negative parameter $`b<0`$ taking the values: (a) $`b=0.40`$; (b) $`b=0.49`$; (c) $`b=0.6`$; (d) $`b=1`$.

Fig. 3. The population difference $`s(t)`$ (dashed line) and phase difference $`x(t)`$ (solid line) as functions of dimensionless time for $`\epsilon =0`$ and the initial conditions $`s_0=0.8`$, $`x_0=0`$, with varying the parameter $`b`$ as: (a) $`b=1`$; (b) $`b=0.85`$; (c) $`b=0.8001`$ (slightly below the critical point $`b_{c1}`$); (d) $`b=0.7999`$ (slightly above the critical point $`b_{c1}`$); (e) $`b=0.5`$; (f) $`b=0.18`$; (g) $`b=0.1999`$ (slightly below the critical point $`b_{c2}`$); (h) $`b=0.2001`$ (slightly above the critical point $`b_{c2}`$); (i) $`b=0.5`$; (j) $`b=1`$.

Fig. 4. Temporal behaviour of $`s(t)`$ (dashed line) and $`x(t)`$ (solid line) for the initial conditions $`s_0=1`$, $`x_0=0`$, fixed $`b=1`$, and a nonzero detuning: (a) $`\epsilon =0.05`$; (b) $`\epsilon =0.1`$.

Fig. 5. Phase portrait on the plane $`\{s,x\}`$ for a finite detuning $`\epsilon =0.1`$ and different values of the pumping-field amplitude: (a) $`b=0.51`$; (b) $`b=0.8`$.


format: text

text:
# Inhomogeneous CuO6 Tilt Distribution and Charge/Spin Correlations in La2-x-yNdySrxCuO4 around commensurate hole concentration

\[

## Abstract

Phononic and magnetic Raman scattering are studied in La2-x-yNdySrxCuO4 with three doping concentrations: $`x1/8`$, $`y=0`$; $`x1/8`$, $`y=0.4`$; and $`x=0.01`$, $`y=0`$. We observe strong disorder in the tilt pattern of the CuO6 octahedra in both the orthorhombic and tetragonal phases which persist down to 10 K and are coupled to bond disorder in the cation layers around $`1/8`$ doping independent of Nd concentration. The weak magnitude of existing charge/spin modulations in the Nd doped structure does not allow us to detect the specific Raman signatures on lattice dynamics or two-magnon scattering around 2200 cm-1.

\]

The evolution of ground state properties with carrier doping and the nature of electronic excitations in correlated systems continues to be a subject of intensive research . After emerging as the mean field solution for the ground state at low hole concentrations in doped Mott-Hubbard antiferromagnets (AF) , stripes have become one paradigm for understanding the properties of high-Tc cuprates . Extensive studies of incommensurate neutron scattering peaks have been performed on La2-xSrxCuO4 (LSCO) crystals and the development of magnetic fluctuations were related to the onset of superconducting order . Transport measurements demonstrate that the charge degrees of freedom self-organize into quasi-1D structures in these materials . It was shown that Nd substitution in LSCO and the proximity to the commensurate $`1/8`$ doping level enhances the stripe correlations and static stripes were observed in La1.48Nd0.4Sr0.12CuO4 (LNSCO) . X-ray diffraction data also observed weak charge order superlattice peaks in this compound .

Lanthanum cuprates are also characterized by disorder introduced by cation substitution and distortions of the CuO planes . Above room temperature LNSCO undergoes a transition from a high temperature tetragonal (HTT) to a low temperature orthorhombic (LTO) phase and around 70 K a transition from the LTO to a low temperature tetragonal (LTT) phase . The LSCO system does not enter the LTT phase but the orthorhombicity decreases with Sr doping. The order parameter of these transitions is the tilt angle of the CuO6 octahedra shown in Fig. 1. It was noticed in LNSCO that the LTO-LTT transition takes place over a range of temperatures and that disorder in the striped phase leads to a glassy nature of the ground state . The coexistence in L(N)SCO of several phases in a complex mixture was suggested by transmission electron microscopy .

The importance of the local CuO6 tilt dynamics for the superconducting properties characterized in cuprates by a short coherence length cannot be ignored. The stabilization of the LTT phase was observed to trace the suppression of superconductivity in Nd doped LSCO  and also in the related La2-xBaxCuO4 compound . A critical value of the tilt angle was associated with the stabilization of magnetic against superconducting order . Rapid suppression of superconductivity, similar to that due to Cu substitution by non-magnetic impurities, was observed with increasing the cation radius variance . By providing information about the energy, scattering width, intensity as well as symmetry of the excitations, Raman scattering has been powerful in studying the structural, electronic, and magnetic properties of cuprates . Optical phonons can be used as local probes of fast changes in the charge distribution and magnetic Raman scattering provides information about local AF correlations . This was the case in related La2-xSrxNiO4+δ compounds where the existence of fully developed striped phases is well established .

Our study provides direct spectroscopic information about the LTO-LTT transition in LNSCO and local deviations from the average structure existent in Nd doped and Nd free LSCO structures. Strong local disorder in the octahedra tilts persists down to T = 10 K and reveals strong disorder in the cation-oxygen layers. The distinct Raman signatures accompanying a transition to a state with deep spin/charge modulations are not observed in the temperature dependence of the two-magnon (2M) scattering around 2200 cm-1 and the $`c`$-axis polarized phonons below 500 cm-1.

Single crystals of LNSCO, $`x=0.01`$ and $`x1/8`$ LSCO were grown by the traveling-solvent floating-zone technique. Raman spectra were taken from the $`(ac)`$ and $`(ab)`$ faces of $`x1/8`$ L(N)SCO and from the $`(a^{}c)`$ face of $`x=0.01`$ LSCO crystals as determined by X-ray diffraction. Unless otherwise stated the laser excitation energy was $`\omega _{in}=1.92`$ eV.

Raman spectra taken in the $`(ca)`$ geometry provide direct information about tetragonal to orthorhombic distortions. In this polarization we probe phononic modes with B2g and B3g symmetries in the LTO phase which become degenerate with Eg symmetry in the LTT phase. The temperature dependence of the modes around 250 cm-1 corresponding to the apical oxygen vibrations parallel to the Cu-O plane in LNSCO is shown in Fig. 2. We observe a broad peak around 245 cm-1 at room temperature which, with cooling, becomes resolved into two components. A new central peak can be seen at 50 K around 248 cm-1 which gains spectral weight as the temperature is decreased to 10 K. While the total integrated intensity of the modes remains constant we observe a redistribution of spectral weight among the three modes as a function of temperature. The split components become weaker but can still be seen as ’orthorhombic satellites’ of the central peak down to 10 K. The coalescence of the features into the 248 cm-1 mode signals the recurrence of a phase with tetragonal symmetry which should be the expected LTT phase of LNSCO. However, the finite residual intensity of the satellites appearing on the tail of the broad central peak shows an incompletely developed LTT phase and that even at 10 K there exists about 7% LTO phase (determined from the relative ratio of phonon intensities).

Raman data in $`(cc)`$ polarization is well suited for the study of lattice dynamics due to weaker coupling to underlying electronic excitations. Group theory predicts five fully symmetric modes at $`𝐤=0`$ in each of the LTO and LTT phases and this is what we observe in Fig. 3 for L(N)SCO crystals in $`(cc)`$ geometry. These modes correspond to: the tilt of the CuO6 octahedra (mode A), the vibration of La/Sr/Nd atoms in the direction imposed by the CuO6 tilt (mode B) and along the $`c`$ axis (mode C), the $`c`$ axis vibrations of the in-plane O (mode D) and apical O atoms (mode E)  (the modes C and E are the two fully symmetric modes predicted by group theory for the HTT phase). The above qualitative description indicates that we could expect a strong coupling between the lowest energy modes (A and B). These two modes can be distinguished because they remain much broader even at the lowest temperature in comparison with the the modes C, D and E which harden and sharpen smoothly through the LTO-LTT transition. The temperature variation of modes C and E inherited from the HTT phase is not as pronounced. The large variation in the energy and width of mode A above the transition, the softening below 70 K, as well as its energy around 110 cm-1 in agreement with neutron scattering studies  show that this mode corresponding to the octahedra tilt is the soft mode of the structural changes. The apparently smooth decrease in energy in the LTT phase is peculiar because this space group is not a subgroup of the LTO group and as a result a true LTO-LTT transition is expected to be of first order. Although unresolved due to broadening effects, the large width of mode A around 70 K shows the coexistence of the LTO and LTT tilts, the latter appearing as a result of folding of the LTO $`Z`$-point to the $`\mathrm{\Gamma }`$-point of the LTT phase as observed for instance in La2NiO4 .

To understand the surprising behavior of the tilt pattern as reflected in the LNSCO phononic data a comparison with different materials from the same class is necessary. In the lower panel of Fig. 3 we show the 10 K $`(cc)`$ polarized Raman spectrum of L(N)SCO. In particular, mode A at 128 cm-1 for $`x=0.01`$ LSCO has a full width at half maximum (FWHM) of 2.5 cm-1. Note in the inset the strongly temperature dependent intensity and width which is a characteristic of a soft mode. For $`x1/8`$ LSCO the same phonon is around 85 cm-1 and its FWHM of about 23.5 cm-1 is larger than the width of mode A in the Nd doped crystal where it is slightly below 20 cm-1. Comparison of the relative phononic widths for both the soft mode and the $`c`$-axis vibrations shows that Nd doping of LSCO crystals and the closer proximity to the T phase induced by Nd doping in the La2CuO4 structure cannot be responsible for the large observed broadening effects. Intrinsic phonon anharmonicity would lead to a broad mode A in $`x=0.01`$ LSCO which is not the case. The tilt disorder across twin domains cannot be the cause of such dramatic effects because the volume fraction occupied by these boundaries is expected to be very small . The 7% relative ratio of the orthorhombic satellites to the central peak in Fig. 2 may be consistent with such a small contribution. If the satellites are indeed due to twinning effects the data show that at 10 K the larger LTT domains are separated by regions of pure LTO tilt. The absence of the broadening effects on the vibrations along the $`c`$-axis points towards an ’anisotropic’ disorder relating primarily to bond randomness along directions parallel to the CuO planes.

Could the spin-lattice coupling or the interaction with the stripe-ordered carriers in CuO planes be the main cause of broadening? Stripe correlations are enhanced in LNSCO which displays however a smaller width of mode A. Also, it is not clear how the $`c`$-axis polarizability would couple to quasi-2D magnetic properties of the CuO planes and why only modes A and B are affected by this interaction. In this sense one expects the in-plane atomic movements to be more sensitive to stripe ordering but we see no similar effects on mode D. Although less probable, spin-lattice induced broadening cannot be completely ruled out and the answer to this question lies in a Sr doping dependence of the $`(cc)`$ polarized spectra. Our data can be reconciled however with recent studies of local structure in Nd free and Nd doped L(N)SCO system . Model analysis of the pair distribution function from X-ray absorption fine structure suggest that in this material class the average structure determined by diffraction is different than the local pattern which is characterized by disorder in the CuO6 tilt direction and magnitude . The Raman data shown in Fig. 3 is a direct spectroscopic evidence that the L(N)SCO system is characterized by strong disorder in the cation layers. The locally fluctuating octahedra tilt is responsible for the observed effects.

Information about the relative magnitude of charge disproportionation in LNSCO can be gained by comparison with Raman spectra in compounds where charge modulations are well established. New Raman active modes have been observed in $`x=0.33`$ and 0.225 La2-xSrxNiO4 by Raman scattering . Lowering of the crystal symmetry at the stripe ordering transition gives rise to folding of the Brillouin zone (BZ) and the appearance of extra $`𝐤=0`$ phononic modes. Charge localization creates non-equivalent Ni sites generating phonon ’splitting’. The $`c`$ axis stretching modes corresponding to La and apical oxygens are split by 14 and 30 cm-1 respectively . Within the 3 cm-1 resolution imposed by the phononic widths we do not observe such a split in our spectra. The ratio of the integrated intensities of the split oxygen modes in Ref.  is about the same as the ratio of doped versus undoped Ni sites. If we assume the same relation holds for the case of cuprates, a factor of $`1/8`$ in split phononic intensity should have been seen in the spectra . We conclude that any charge ordering taking place in our case is much weaker than in the related compounds referred to above. This does not contradict the X-ray diffraction data  which estimated a factor of 102 relative difference between the magnitude of charge modulations in cuprates and nickelates.

For 2D magnetic square lattices a 2M peak in the $`(a^{}b^{})`$ geometry is expected . Fig. 4 shows 2M scattering around 2200 cm-1 at 300 and 10 K taken with the resonant $`\omega _{in}`$ = 3.05 eV incident frequency. As in other tetragonal 2D AFs we observe the spin pair excitations in the expected geometry which corresponds to B1g symmetry in the HTT phase. Renormalization of the 2M excitation energy occurs across the spin ordering transition in the striped phase of La2-xSrxNiO4. The 2M Raman band around 1650 cm-1 characteristic of the undoped case  is not seen while the spectra exhibit lower energy peaks below the magnetic ordering temperature as a result of the new spin exchange channels within and across the antiphase AF domains . Our data in LNSCO shows only slight changes from 300 to 10 K emphasizing weak local spin modulations.

The differences we observe between cuprates and nickelates can be related to the much stronger carrier self-confinement in the latter . It has also been shown  that anomalies in phonon dispersions occur in La2-xSrxCuO4 at points in the BZ commensurate with charge ordering wavevectors inferred from neutron scattering studies. But as discussed, the charge modulation in Nd doped structures, where the stripe correlations were shown to be stabilized, is too weak to produce observable changes (at $`𝐤=0`$) in the lattice unit cell on time scales longer than the phononic frequencies. The number of phononic modes we observe can be explained solely in terms of LTO/LTT distortions. Our data, however, does not contradict the possible existence of charge modulations in the Cu-O plane. In fact, the dynamics in the cation-O layers and the magnitude of octahedra tilt disorder affects the carrier distribution and our Raman results impose constraints on the magnitude of the charge modulations.

We studied lattice dynamics and magnetic scattering in stripe ordered L(N)SCO crystals. The transition to the LTT phase in LNSCO is not complete as shown by the behavior of the split oxygen modes around 250 cm-1. The large widths of the soft modes demonstrate a locally fluctuating pattern of the CuO6 octahedra tilts. This disorder effect is even stronger in the Nd free LSCO compound at $`x1/8`$ while it is not present for $`x=0.01`$ Sr concentration. Quantitative comparison with lanthanum-nickelate compounds shows that the magnitude of the existent spin/charge disproportionations in L(N)SCO are relatively weak and that, as in other cuprates, the two-magnon excitation can be described within the usual high temperature tetragonal structure.


text:
# Finding and evaluating community structure in networks

## I Introduction

Empirical studies and theoretical modeling of networks have been the subject of a large body of recent research in statistical physics and applied mathematics Strogatz01 ; AB02 ; DM03b ; Newman03d . Network ideas have been applied with great success to topics as diverse as the Internet and the world wide web FFF99 ; AJB99 ; Broder00 , epidemiology KG99 ; MN00a ; PV01a ; ML01 , scientific citation and collaboration Redner98 ; Newman01a , metabolism Jeong00 ; WF01 , and ecosystems DWM02a ; CGA02 , to name but a few. A property that seems to be common to many networks is community structure, the division of network nodes into groups within which the network connections are dense, but between which they are sparser—see Fig. 1. The ability to find and analyze such groups can provide invaluable help in understanding and visualizing the structure of networks. In this paper we show how this can be achieved.

The study of community structure in networks has a long history. It is closely related to the ideas of graph partitioning in graph theory and computer science, and hierarchical clustering in sociology GJ79 ; Scott00 . Before presenting our own findings, it is worth reviewing some of this preceding work, to understand its achievements and where it falls short.

Graph partitioning is a problem that arises in, for example, parallel computing. Suppose we have a number $`n`$ of intercommunicating computer processes, which we wish to distribute over a number $`g`$ of computer processors. Processes do not necessarily need to communicate with all others, and the pattern of required communications can be represented by a graph or network in which the vertices represent processes and edges join process pairs that need to communicate. The problem is to allocate the processes to processors in such a way as roughly to balance the load on each processor, while at the same time minimizing the number of edges that run between processors, so that the amount of interprocessor communication (which is normally slow) is minimized. In general, finding an exact solution to a partitioning task of this kind is believed to be an NP-complete problem, making it prohibitively difficult to solve for large graphs, but a wide variety of heuristic algorithms have been developed that give acceptably good solutions in many cases, the best known being perhaps the Kernighan–Lin algorithm KL70 , which runs in time $`\mathrm{O}(n^3)`$ on sparse graphs.

A solution to the graph partitioning problem is however not particularly helpful for analyzing and understanding networks in general. If we merely want to find if and how a given network breaks down into communities, we probably don’t know how many such communities there are going to be, and there is no reason why they should be roughly the same size. Furthermore, the number of inter-community edges needn’t be strictly minimized either, since more such edges are admissible between large communities than between small ones.

As far as our goals in this paper are concerned, a more useful approach is that taken by social network analysis with the set of techniques known as hierarchical clustering. These techniques are aimed at discovering natural divisions of (social) networks into groups, based on various metrics of similarity or strength of connection between vertices. They fall into two broad classes, agglomerative and divisive Scott00 , depending on whether they focus on the addition or removal of edges to or from the network. In an agglomerative method, similarities are calculated by one method or another between vertex pairs, and edges are then added to an initially empty network ($`n`$ vertices with no edges) starting with the vertex pairs with highest similarity. The procedure can be halted at any point, and the resulting components in the network are taken to be the communities. Alternatively, the entire progression of the algorithm from empty graph to complete graph can be represented in the form of a tree or dendrogram such as that shown in Fig. 2. Horizontal cuts through the tree represent the communities appropriate to different halting points.

Agglomerative methods based on a wide variety of similarity measures have been applied to different networks. Some networks have natural similarity metrics built in. For example, in the widely studied network of collaborations between film actors WS98 ; ASBS00 , in which two actors are connected if they have appeared in the same film, one could quantify similarity by how many films actors have appeared in together ML00 . Other networks have no natural metric, but suitable ones can be devised using correlation coefficients, path lengths, or matrix methods. A well known example of an agglomerative clustering method is the Concor algorithm of Breiger et al. BBA75 .

Agglomerative methods have their problems however. One concern is that they fail with some frequency to find the correct communities in networks were the community structure is known, which makes it difficult to place much trust in them in other cases. Another is their tendency to find only the cores of communities and leave out the periphery. The core nodes in a community often have strong similarity, and hence are connected early in the agglomerative process, but peripheral nodes that have no strong similarity to others tend to get neglected, leading to structures like that shown in Fig. 3. In this figure, there are a number of peripheral nodes whose community membership is obvious to the eye—in most cases they have only a single link to a specific community—but agglomerative methods often fail to place such nodes correctly.

In this paper, therefore, we focus on divisive methods. These methods have been relatively little studied in the previous literature, either in social network theory or elsewhere, but, as we will see, seem to offer a lot of promise. In a divisive method, we start with the network of interest and attempt to find the *least* similar connected pairs of vertices and then remove the edges between them. By doing this repeatedly, we divide the network into smaller and smaller components, and again we can stop the process at any stage and take the components at that stage to be the network communities. Again, the process can be represented as a dendrogram depicting the successive splits of the network into smaller and smaller groups.

The approach we take follows roughly these lines, but adopts a somewhat different philosophical viewpoint. Rather than looking for the most weakly connected vertex pairs, our approach will be to look for the edges in the network that are most “between” other vertices, meaning that the edge is, in some sense, responsible for connecting many pairs of others. Such edges need not be weak at all in the similarity sense. How this idea works out in practice will become clear in the course of the presentation.

Briefly then, the outline of this paper is as follows. In Sec. II we describe the crucial concepts behind our methods for finding community structure in networks and show how these concepts can be turned into a concrete prescription for performing calculations. In Sec. III we describe in detail the implementation of our methods. In Sec. IV we consider ways of determining when a particular division of a network into communities is a good one, allowing us to quantify the success of our community-finding algorithms. And in Sec. V we give a number of applications of our algorithms to particular networks, both real and artificial. In Sec. VI we give our conclusions. A brief report of some of the work contained in this paper has appeared previously as Ref. GN02, .

## II Finding communities in a network

In this paper we present a class of new algorithms for network clustering, i.e., the discovery of community structure in networks. Our discussion focuses primarily on networks with only a single type of vertex and a single type of undirected, unweighted edge, although generalizations to more complicated network types are certainly possible.

There are two central features that distinguish our algorithms from those that have preceded them. First, our algorithms are divisive rather than agglomerative. Divisive algorithms have occasionally been studied in the past, but, as discussed in the introduction, ours differ in focusing not on removing the edges between vertex pairs with lowest similarity, but on finding edges with the highest “betweenness,” where betweenness is some measure that favors edges that lie between communities and disfavors those that lie inside communities.

To make things more concrete, we give some examples of the types of betweenness measures we will be looking at. All of them are based on the same idea. If two communities are joined by only a few inter-community edges, then all paths through the network from vertices in one community to vertices in the other must pass along one of those few edges. Given a suitable set of paths, one can count how many go along each edge in the graph, and this number we then expect to be largest for the inter-community edges, thus providing a method for identifying them. Our different measures correspond to various implementations of this idea.

1. The simplest example of such a betweenness measure is that based on shortest (geodesic) paths: we find the shortest paths between all pairs of vertices and count how many run along each edge. To the best of our knowledge this measure was first introduced by Anthonisse in a never-published technical report in 1971 Anthonisse71 . Anthonisse called it “rush,” but we prefer the term edge betweenness, since the quantity is a natural generalization to edges of the well-known (vertex) betweenness measure of Freeman Freeman77 , which was the inspiration for our approach. When we need to distinguish it from the other betweenness measures considered in this paper, we will refer to it as shortest-path betweenness. A fast algorithm for calculating the shortest-path betweenness is given in Sec. III.1.
2. The shortest-path betweenness can be thought of in terms of signals traveling through a network. If signals travel from source to destination along geodesic network paths, and all vertices send signals at the same constant rate to all others, then the betweenness is a measure of the rate at which signals pass along each edge. Suppose however that signals do not travel along geodesic paths, but instead just perform a random walk about the network until they reach their destination. This gives us another measure on edges, the random-walk betweenness: we calculate the expected net number of times that a random walk between a particular pair of vertices will pass down a particular edge and sum over all vertex pairs. The random-walk betweenness can be calculated using matrix methods, as described in Sec. III.3.
3. Another betweenness measure is motivated by ideas from elementary circuit theory. We consider the circuit created by placing a unit resistance on each edge of the network and unit current source and sink at a particular pair of vertices. The resulting current flow in the network will travel from source to sink along a multitude of paths, those with least resistance carrying the greatest fraction of the current. The current-flow betweenness for an edge we define to be the absolute value of the current along the edge summed over all source/sink pairs. It can be calculated using Kirchhoff’s laws, as described in Sec. III.2. In fact, as we will show, the current-flow betweenness turns out to be exactly the same as the random walk betweenness of the previous paragraph, but we nonetheless consider it separately since it leads to a simpler derivation of the measure.

These measures are only suggestions; many others are possible and may well be appropriate for specific applications. Measures (1) and (2) are in some sense extremes in the spectrum of possibilities, one corresponding to signals that know exactly where they are going, and the other to signals that have no idea where they are going. As we will see, however, these two measures actually give rather similar results, indicating that the precise choice of betweenness measure may not, at least for the types of applications considered here, be that important.

The second way in which our methods differ from previous ones is in the inclusion of a “recalculation step” in the algorithm. If we were to perform a standard divisive clustering based on edge betweenness we would calculate the edge betweenness for all edges in the network and then remove edges in decreasing order of betweenness to produce a dendrogram like that of Fig. 2, showing the order in which the network split up.

However, once the first edge in the network is removed in such an algorithm, the betweenness values for the remaining edges will no longer reflect the network as it now is. This can give rise to unwanted behaviors. For example, if two communities are joined by two edges, but, for one reason or another, most paths between the two flow along just one of those edges, then that edge will have a high betweenness score and the other will not. An algorithm that calculated betweennesses only once and then removed edges in betweenness order would remove the first edge early in the course of its operation, but the second might not get removed until much later. Thus the obvious division of the network into two parts might not be discovered by the algorithm. In the worst case the two parts themselves might be individually broken up before the division between the two is made. In practice, problems like this crop up in real networks with some regularity and render algorithms of this type ineffective for the discovery of community structure.

The solution, luckily, is obvious. We simply recalculate our betweenness measure after the removal of each edge. This certainly adds to the computational effort of performing the calculation, but its effect on the results is so desirable that we consider the price worth paying.

Thus the general form of our community structure finding algorithm is as follows:

1. Calculate betweenness scores for all edges in the network.
2. Find the edge with the highest score and remove it from the network.
3. Recalculate betweenness for all remaining edges.
4. Repeat from step 2.

In fact, it appears that the recalculation step is the most important feature of the algorithm, as far as getting satisfactory results is concerned. As mentioned above, our studies indicate that, once one hits on the idea of using betweenness measures to weight edges, the exact measure one uses appears not to influence the results highly. The recalculation step, on the other hand, is absolutely crucial to the operation of our methods. This step was missing from previous attempts at solving the clustering problem using divisive algorithms, and yet without it the results are very poor indeed, failing to find known community structure even in the simplest of cases. In Sec. V.2 we give an example comparing the performance of the algorithm on a particular network with and without the recalculation step.

In the following sections we discuss implementation and give examples of our algorithms for finding community structure. For the reader who merely wants to know what algorithm they should use for their own problem, let us give an immediate answer: for most problems, we recommend the algorithm with betweenness scores calculated using the shortest-path betweenness measure (1) above. This measure appears to work well and is the quickest to calculate—as described in Sec. III.1, it can be calculated for all edges in time $`\mathrm{O}(mn)`$, where $`m`$ is the number of edges in the graph and $`n`$ is the number of vertices. This is the only version of the algorithm that we discussed in Ref. GN02,  111Following the publication of that paper, the algorithm has been implemented in the commercial network analysis software package UCInet and in the open-source network library JUNG. (See http://www.analytictech.com/ and http://jung.sourceforge.net/ respectively.). The other versions we discuss, while being of some pedagogical interest, make greater computational demands, and in practice seem to give results no better than the shortest-path method.

## III Implementation

In theory, the descriptions of the last section completely define the methods we consider in this paper, but in practice there are a number of tricks to their implementation that are important for turning the description into a workable computer algorithm.

Essentially all of the work in the algorithm is in the calculation of the betweenness scores for the edges; the job of finding and removing the highest-scoring edge is trivial and not computationally demanding. Let us tackle our three suggested betweenness measures in turn.

### III.1 Shortest-path betweenness

At first sight, it appears that calculating the edge betweenness measure based on geodesic paths for all edges will take $`\mathrm{O}(mn^2)`$ operations on a graph with $`m`$ edges and $`n`$ vertices: calculating the shortest path between a particular pair of vertices can be done using breadth-first search in time $`\mathrm{O}(m)`$ AMO93 ; CLRS01 , and there are $`\mathrm{O}(n^2)`$ vertex pairs. Recently however new algorithms have been proposed by Newman Newman01c  and independently by Brandes Brandes01  that can perform the calculation faster than this, finding all betweennesses in $`\mathrm{O}(mn)`$ time. Both Newman and Brandes gave algorithms for the standard Freeman vertex betweenness, but it is trivial to adapt their algorithms for edge betweenness. We describe the resulting method here for the algorithm of Newman.

Breadth-first search can find shortest paths from a single vertex $`s`$ to all others in time $`\mathrm{O}(m)`$. In the simplest case, when there is only a single shortest path from the source vertex to any other (we will consider other cases in a moment) the resulting set of paths forms a shortest-path tree—see Fig. 4a. We can now use this tree to calculate the contribution to betweenness for each edge from this set of paths as follows. We find first the “leaves” of the tree, i.e., those nodes such that no shortest paths to other nodes pass through them, and we assign a score of 1 to the single edge that connects each to the rest of the tree, as shown in the figure. Then, starting with those edges that are farthest from the source vertex on the tree, i.e., lowest in Fig. 4a, we work upwards, assigning a score to each edge that is 1 plus the sum of the scores on the neighboring edges immediately below it. When we have gone though all edges in the tree, the resulting scores are the betweenness counts for the paths from vertex $`s`$. Repeating the process for all possible vertices $`s`$ and summing the scores, we arrive at the full betweenness scores for shortest paths between all pairs. The breadth-first search and the process of working up through the tree both take worst-case time $`\mathrm{O}(m)`$ and there are $`n`$ vertices total, so the entire calculation takes time $`\mathrm{O}(mn)`$ as claimed.

This simple case serves to illustrate the basic principle behind the algorithm. In general, however, it is not the case that there is only a single shortest path between any pair of vertices. Most networks have at least some vertex pairs between which there are several geodesic paths of equal length. Figure 4b shows a simple example of a shortest path “tree” for a network with this property. The resulting structure is in fact no longer a tree, and in such cases an extra step is required in the algorithm to calculate the betweenness correctly.

In the traditional definition of vertex betweenness Freeman77  multiple shortest paths between a pair of vertices are given equal weights summing to 1. For example, if there are three shortest paths, each will be given weight $`\frac{1}{3}`$. We adopt the same definition for our edge betweenness (as did Anthonisse in his original work Anthonisse71 , although other definitions are possible GKK01b ). Note that the paths may run along the same edge or edges for some part of their length, resulting in edges with greater weight. To calculate correctly what fraction of the paths flow along each edge in the network, we generalize the breadth-first search part of the calculation, as follows.

Consider Fig. 4b and suppose we are performing a breadth-first search starting at vertex $`s`$. We carry out the following steps:

1. The initial vertex $`s`$ is given distance $`d_s=0`$ and a weight $`w_s=1`$.
2. Every vertex $`i`$ adjacent to $`s`$ is given distance $`d_i=d_s+1=1`$, and weight $`w_i=w_s=1`$.
3. For each vertex $`j`$ adjacent to one of *those* vertices $`i`$ we do one of three things:

1. If $`j`$ has not yet been assigned a distance, it is assigned distance $`d_j=d_i+1`$ and weight $`w_j=w_i`$.
    2. If $`j`$ has already been assigned a distance and $`d_j=d_i+1`$, then the vertex’s weight is increased by $`w_i`$, that is $`w_jw_j+w_i`$.
    3. If $`j`$ has already been assigned a distance and $`d_j<d_i+1`$, we do nothing.
4. Repeat from step 3 until no vertices remain that have assigned distances but whose neighbors do not have assigned distances.

In practice, this algorithm can be implemented most efficiently using a queue or first-in/first-out buffer to store the vertices that have been assigned a distance, just as in the standard breadth-first search.

Physically, the weight on a vertex $`i`$ represents the number of distinct paths from the source vertex to $`i`$. These weights are precisely what we need to calculate our edge betweennesses, because if two vertices $`i`$ and $`j`$ are connected, with $`j`$ farther than $`i`$ from the source $`s`$, then the fraction of a geodesic path from $`j`$ through $`i`$ to $`s`$ is given by $`w_i/w_j`$. Thus, to calculate the contribution to edge betweenness from all shortest paths starting at $`s`$, we need only carry out the following steps:

1. Find every “leaf” vertex $`t`$, i.e., a vertex such that no paths from $`s`$ to other vertices go though $`t`$.
2. For each vertex $`i`$ neighboring $`t`$ assign a score to the edge from $`t`$ to $`i`$ of $`w_i/w_t`$.
3. Now, starting with the edges that are farthest from the source vertex $`s`$—lower down in a diagram such as Fig. 4b—work up towards $`s`$. To the edge from vertex $`i`$ to vertex $`j`$, with $`j`$ being farther from $`s`$ than $`i`$, assign a score that is 1 plus the sum of the scores on the neighboring edges immediately below it (i.e., those with which it shares a common vertex), all multiplied by $`w_i/w_j`$.
4. Repeat from step 3 until vertex $`s`$ is reached.

Now repeating this process for all $`n`$ source vertices $`s`$ and summing the resulting scores on the edges gives us the total betweenness for all edges in time $`\mathrm{O}(mn)`$.

We now have to repeat this calculation for each edge removed from the network, of which there are $`m`$, and hence the complete community structure algorithm based on shortest-path betweenness operates in worst-case time $`\mathrm{O}(m^2n)`$, or $`\mathrm{O}(n^3)`$ time on a sparse graph. In our experience, this typically makes it tractable for networks of up to about $`n=\mathrm{10\hspace{0.17em}000}`$ vertices, with current (circa 2003) desktop computers. In some special cases one can do better. In particular, we note that the removal of an edge only affects the betweenness of other edges that fall in the same component, and hence that we need only recalculate betweennesses in that component. Networks with strong community structure often break apart into separate components quite early in the progress of the algorithm, substantially reducing the amount of work that needs to be done on subsequent steps. Whether this results in a change in the computational complexity of the algorithm for any commonly occurring classes of graphs is an open question, but it certainly gives a substantial speed boost to many of the calculations described in this paper.

Some networks are directed, i.e., their edges run in one direction only. The world wide web is an example; links in the web point in one direction only from one web page to another. One could imagine a generalization of the shortest-path betweenness that allowed for directed edges by counting only those paths that travel in the forward direction along edges. Such a calculation is a trivial variation on the one described above. However, we have found that in many cases it is better to ignore the directed nature of a network in calculating community structure. Often an edge acts simply as an indication of a connection between two nodes, and its direction is unimportant. For example, in Ref. GN02,  we applied our algorithm to a food web of predator-prey interactions between marine species. Predator-prey interactions are clearly directed—one species may eat another, but it is unlikely that the reverse is simultaneously true. However, as far as community structure is concerned, we want to know only which species have interactions with which others. We find, therefore, that our algorithm applied to the undirected version of the food web works well at picking out the community structure, and no special algorithm is needed for the directed case. We give another example of our method applied to a directed graph in Sec. V.4.

### III.2 Resistor networks

As examples of betweenness measures that take more than just shortest paths into account, we proposed in Sec. II measures based on random walks and resistor networks. In fact, as we now show, when appropriately defined these two measures are precisely the same. Here we derive the resistance measure first, since it turns out to be simpler; in the following section we derive the random walk measure and show that the two are equivalent.

Consider the network created by placing a unit resistance on every edge of our network, a unit current source at vertex $`s`$, and a unit current sink at vertex $`t`$ (see Fig. 5). Clearly the current between $`s`$ and $`t`$ will flow primarily along short paths, but some will flow along longer ones, roughly in inverse proportion to their length. We will use the absolute magnitude of the current flow as our betweenness score for each source/sink pair.

The current flows in the network are governed by Kirchhoff’s laws. To solve them we proceed as follows for each separate component of the graph. Let $`V_i`$ be the voltage at vertex $`i`$, measured relative to any convenient point. Then for all $`i`$ we have

$$\underset{j}{}A_{ij}(V_iV_j)=\delta _{is}\delta _{it},$$
(1)
where $`A_{ij}`$ is the $`ij`$ element of the adjacency matrix of the graph, i.e., $`A_{ij}=1`$ if $`i`$ and $`j`$ are connected by an edge and $`A_{ij}=0`$ otherwise. The left-hand side of Eq. (1) represents the net current flow out of vertex $`i`$ along edges of the network, and the right-hand side represents the source and sink. Defining $`k_i=_jA_{ij}`$, which is the vertex degree, and creating a diagonal matrix $`𝐃`$ with these degrees on the diagonal $`D_{ii}=k_i`$, this equation can be written in matrix form as $`(𝐃𝐀)𝐕=𝐬`$, where the source vector $`𝐬`$ has components

$$s_i=\{\begin{array}{cc}\hfill +1& \text{for }i=s\hfill \\ \hfill 1& \text{for }i=t\hfill \\ \hfill 0& \text{otherwise.}\hfill \end{array}$$
(2)
We cannot directly invert the matrix $`𝐃𝐀`$ to get the voltage vector $`𝐕`$, because the matrix (which is just the graph Laplacian) is singular. This is equivalent to saying that there is one undetermined degree of freedom corresponding to the choice of reference potential for measuring the voltages. We can add any constant to a solution for the vertex voltages and get another solution—only the voltage differences matter. In choosing the reference potential, we fix this degree of freedom, leaving only $`n1`$ more to be determined. In mathematical terms, once any $`n1`$ of the equations in our matrix formulation are satisfied, the remaining one is also automatically satisfied so long as current is conserved in the network as a whole, i.e., so long as $`_is_i=0`$, which is clearly true in this case.

Choosing any vertex $`v`$ to be the reference point, therefore, we remove the row and column corresponding to that vertex from $`𝐃`$ and $`𝐀`$ before inverting. Denoting the resulting $`(n1)\times (n1)`$ matrices $`𝐃_v`$ and $`𝐀_v`$, we can then write

$$𝐕=(𝐃_v𝐀_v)^1𝐬.$$
(3)
Calculation of the currents in the network thus involves inverting $`𝐃_v𝐀_v`$ once for any convenient choice of $`v`$, and taking the differences of pairs of columns to get the voltage vector $`𝐕`$ for each possible source/sink pair. (The voltage for the one missing vertex $`v`$ is always zero, by hypothesis.) The absolute magnitudes of the differences of voltages along each edge give us betweenness scores for the given source and sink. Summing over all sources and sinks, we then get our complete betweenness score.

The matrix inversion takes time $`\mathrm{O}(n^3)`$ in the worst case, while the subsequent calculation of betweennesses takes time $`\mathrm{O}(mn^2)`$, where as before $`m`$ is the number of edges and $`n`$ the number of vertices in the graph. Thus, the entire community structure algorithm, including the recalculation step, will take $`\mathrm{O}\left((n+m)mn^2\right)`$ time to complete, or $`\mathrm{O}(n^4)`$ on a sparse graph. Although, as we will see, the algorithm is good at finding community structure, this poor performance makes it practical only for smaller graphs; a few hundreds of vertices is the most that we have been able to do. It is for this reason that we recommend using the shortest-path betweenness algorithm in most cases, which gives results about as good or better with considerably less effort.

### III.3 Random walks

The random-walk betweenness described in Sec. II requires us to calculate how often on average random walks starting at vertex $`s`$ will pass down a particular edge from vertex $`v`$ to vertex $`w`$ (or vice versa) before finding their way to a given target vertex $`t`$. To calculate this quantity we proceed as follows for each separate component of the graph.

As before, let $`A_{ij}`$ be an element of the adjacency matrix such that $`A_{ij}=1`$ if vertices $`i`$ and $`j`$ are connected by an edge and $`A_{ij}=0`$ otherwise. Consider a random walk that on each step decides uniformly between the neighbors of the current vertex $`j`$ and takes a step to one of them. The number of neighbors is just the degree of the vertex $`k_j=_iA_{ij}`$, and the probability for the transition from $`j`$ to $`i`$ is $`A_{ij}/k_j`$, which we can regard as an element of the matrix $`𝐌=𝐀𝐃^1`$, where $`𝐃`$ is the diagonal matrix with $`D_{ii}=k_i`$.

We are interested in walks that terminate when they reach the target $`t`$, so that $`t`$ is an absorbing state. The most convenient way to represent this is just to remove entirely the vertex $`t`$ from the graph, so that no walk ever reaches any other vertex from $`t`$. Thus let $`𝐌_t=𝐀_t𝐃_t^1`$ be the matrix $`𝐌`$ with the $`t`$th row and column removed (and similarly for $`𝐀_t`$ and $`𝐃_t`$).

Now the probability that a walk starts at $`s`$, takes $`n`$ steps, and ends up at some other vertex (not $`t`$), is given by the $`is`$ element of $`𝐌_t^n`$, which we denote $`[𝐌_t^n]_{is}`$. In particular, walks end up at $`v`$ and $`w`$ with probabilities $`[𝐌_t^n]_{vs}`$ and $`[𝐌_t^n]_{ws}`$, and of those a fraction $`1/k_v`$ and $`1/k_w`$ respectively then pass along the edge $`(v,w)`$ in one direction or the other. (Note that they may also have passed along this edge an arbitrary number of times before reaching this point.) Summing over all $`n`$, the mean number of times that a walk of any length traverses the edge from $`v`$ to $`w`$ is $`k_v^1[(𝐈𝐌_t)^1]_{vs}`$, and similarly for walks that go from $`w`$ to $`v`$.

To highlight the similarity with the current-flow betweenness of Sec. III.2, let us denote these two numbers $`V_v`$ and $`V_w`$ respectively. Then we can write

$$𝐕=𝐃^1(𝐈𝐌_t)^1𝐬=(𝐃_t𝐀_t)^1𝐬,$$
(4)
where the source vector $`𝐬`$ is the vector whose components are all 0 except for a single 1 in the position corresponding to the source vertex $`s`$.

Now we define our random-walk betweenness for the edge $`(v,w)`$ to be the absolute value of the *difference* of the two probabilities $`V_v`$ and $`V_w`$, i.e., the net number of times the walk passes along the edge in one direction. This seems a natural definition—it makes little sense to accord an edge high betweenness simply because a walk went back and forth along it many times. It is the difference between the numbers of times the edge is traversed in either direction that matters 222In fact we have tried counting each traversal separately, but this method gives extremely poor results, confirming our intuition that this would not be a good betweenness measure..

But now we see that this method is very similar to the resistor network calculation of Sec. III.2. In that calculation we also evaluated $`(𝐃_t𝐀_t)^1𝐬`$ for a suitable source vector and then took differences of the resulting numbers. The only difference is that in the current-flow calculation we had a sink term in $`𝐬`$ as well as a source. Purely for the purposes of mathematical convenience, we can add such a sink in the present case at the target vertex $`t`$—this makes no difference to the solution for $`𝐕`$ since the $`t`$th row has been removed from the equations anyway. By doing this, however, we turn the equations into precisely the form of the current-flow calculation, and hence it becomes clear that the two measures are numerically identical, although their derivation is quite different. (It also immediately follows that we can remove any row or column and still get the same answer—it doesn’t have to be row and column $`t`$, although physically this choice makes the most sense.)

## IV Quantifying the strength of community structure

As we show in Sec. V, our community structure algorithms do an excellent job of recovering known communities both in artificially generated random networks and in real-world examples. However, in practical situations the algorithms will normally be used on networks for which the communities are not known ahead of time. This raises a new problem: how do we know when the communities found by the algorithm are good ones? Our algorithms always produce *some* division of the network into communities, even in completely random networks that have no meaningful community structure, so it would be useful to have some way of saying how good the structure found is. Furthermore, the algorithms’ output is in the form of a dendrogram which represents an entire nested hierarchy of possible community divisions for the network. We would like to know which of these divisions are the best ones for a given network—where we should cut the dendrogram to get a sensible division of the network.

To answer these questions we now define a measure of the quality of a particular division of a network, which we call the modularity. This measure is based on a previous measure of assortative mixing proposed by Newman Newman03c . Consider a particular division of a network into $`k`$ communities. Let us define a $`k\times k`$ symmetric matrix $`𝐞`$ whose element $`e_{ij}`$ is the fraction of all edges in the network that link vertices in community $`i`$ to vertices in community $`j`$ 333As discussed in Newman03c , it is crucial to make sure each edge is counted only once in the matrix $`e_{ij}`$—the same edge should not appear both above and below the diagonal. Alternatively, an edge linking communities $`i`$ and $`j`$ can be split, half-and-half, between the $`ij`$ and $`ji`$ elements, which has the advantage of making the matrix symmetric. Either way, there are a number of factors of 2 in the calculation that must be watched carefully, lest they escape one’s attention and make mischief.. (Here we consider all edges in the original network—even after edges have been removed by the community structure algorithm our modularity measure is calculated using the full network.)

The trace of this matrix $`Tr𝐞=_ie_{ii}`$ gives the fraction of edges in the network that connect vertices in the same community, and clearly a good division into communities should have a high value of this trace. The trace on its own, however, is not a good indicator of the quality of the division since, for example, placing all vertices in a single community would give the maximal value of $`Tr𝐞=1`$ while giving no information about community structure at all.

So we further define the row (or column) sums $`a_i=_je_{ij}`$, which represent the fraction of edges that connect to vertices in community $`i`$. In a network in which edges fall between vertices without regard for the communities they belong to, we would have $`e_{ij}=a_ia_j`$. Thus we can define a modularity measure by

$$Q=\underset{i}{}\left(e_{ii}a_i^2\right)=Tr𝐞𝐞^2,$$
(5)
where $`𝐱`$ indicates the sum of the elements of the matrix $`𝐱`$. This quantity measures the fraction of the edges in the network that connect vertices of the same type (i.e., within-community edges) minus the expected value of the same quantity in a network with the same community divisions but random connections between the vertices. If the number of within-community edges is no better than random, we will get $`Q=0`$. Values approaching $`Q=1`$, which is the maximum, indicate strong community structure 444In Ref. Newman03c, , the measure was normalized by dividing by its value on a network with perfect mixing, so that we always get 1 for such a network. We find however that doing this in the present case masks some of the useful information to be gained from the value of $`Q`$, and hence that it is better to use the unnormalized measure. In general, this unnormalized measure will not reach a value of 1, even on a perfectly mixed network.. In practice, values for such networks typically fall in the range from about $`0.3`$ to $`0.7`$. Higher values are rare.

The expected error on $`Q`$ can be calculated by treating each edge in the network as an independent measurement of the contributions to the elements of the matrix $`𝐞`$. A simple jackknife procedure works well Efron79 ; Newman03c .

Typically, we will calculate $`Q`$ for each split of a network into communities as we move down the dendrogram, and look for local peaks in its value, which indicate particularly satisfactory splits. Usually we find that there are only one or two such peaks and, as we will show in the next section, in cases where the community structure is known beforehand by some means we find that the positions of these peaks correspond closely to the expected divisions. The height of a peak is a measure of the strength of the community division.

## V Applications

In this section we give a number of applications of our algorithms to particular problems, illustrating their operation, and their use in understanding the structure of complex networks.

### V.1 Tests on computer-generated networks

First, as a controlled test of how well our algorithms perform, we have generated networks with known community structure, to see if the algorithms can recognize and extract this structure.

We have generated a large number of graphs with $`n=128`$ vertices, divided into four communities of 32 vertices each. Edges were placed independently at random between vertex pairs with probability $`p_{\mathrm{in}}`$ for an edge to fall between vertices in the same community and $`p_{\mathrm{out}}`$ to fall between vertices in different communities. The values of $`p_{\mathrm{in}}`$ and $`p_{\mathrm{out}}`$ were chosen to make the expected degree of each vertex equal to 16. In Fig. 6 we show a typical dendrogram from the analysis of such a graph using the shortest-path betweenness version of our algorithm. (In fact, for the sake of clarity, the figure is for a 64-node version of the graph.) Results for the random walk version are similar. At the top of the figure we also show the modularity, Eq. (5), for the same calculation, plotted as a function of position in the dendrogram. That is, the graph is aligned with the dendrogram so that one can read off modularity values for different divisions of the network directly. As we can see, the modularity has a single clear peak at the point where the network breaks into four communities, as we would expect. The peak value is around $`0.5`$, which is typical.

In Fig. 7 we show the fraction of vertices in our computer-generated network sample classified correctly into the four communities by our algorithms, as a function of the mean number $`z_{\mathrm{out}}`$ of edges from each vertex to vertices in other communities. As the figure shows, both the shortest-path and random-walk versions of the algorithm perform excellently, with more than 90% of all vertices classified correctly from $`z_{\mathrm{out}}=0`$ all the way to around $`z_{\mathrm{out}}=6`$. Only for $`z_{\mathrm{out}}6`$ does the classification begin to deteriorate markedly. In other words, our algorithm correctly identifies the community structure in the network almost all the way to the point $`z_{\mathrm{out}}=8`$ at which each vertex has on average the same number of connections to vertices outside its community as it does to those inside.

The shortest-path version of the algorithm does however perform noticeably better than the random-walk version, especially for the more difficult cases where $`z_{\mathrm{out}}`$ is large. Given that the random-walk algorithm is also more computationally demanding, there seems little reason to use it rather than the shortest-path algorithm, and hence, as discussed previously, we recommend the latter for most applications. (To be fair, the random-walk algorithm does slightly out-perform the shortest-path algorithm in the example addressed in the following section, although, being only a single case, it is hard to know whether this is significant.)

### V.2 Zachary’s karate club network

We now turn to applications of our methods to real-world network data. Our first such example is taken from one of the classic studies in social network analysis. Over the course of two years in the early 1970s, Wayne Zachary observed social interactions between the members of a karate club at an American university Zachary77 . He constructed networks of ties between members of the club based on their social interactions both within the club and away from it. By chance, a dispute arose during the course of his study between the club’s administrator and its principal karate teacher over whether to raise club fees, and as a result the club eventually split in two, forming two smaller clubs, centered around the administrator and the teacher.

In Fig. 8 we show a consensus network structure extracted from Zachary’s observations before the split. Feeding this network into our algorithms we find the results shown in Fig. 9. In the left-most two panels we show the dendrograms generated by the shortest-path and random-walk versions of our algorithm, along with the modularity measures for the same. As we see, both algorithms give reasonably high values for the modularity when the network is split into two communities—around 0.4 in each case—indicating that there is a strong natural division at this level. What’s more, the divisions in question correspond almost perfectly to the actual divisions in the club revealed by which group each club member joined after the club split up. (The shapes of the vertices representing the two factions are the same as those of Fig. 8.) Only one vertex, vertex 3, is misclassified by the shortest-path version of the method, and none are misclassified by the random-walk version—the latter gets a perfect score on this test. (On the other hand, the two-community split fails to produce a local maximum in the modularity for the random-walk method, unlike the shortest-path method for which there is a local maximum precisely at this point.)

In the last panel of Fig. 9 we show the dendrogram and modularity for an algorithm based on shortest-path betweenness but without the crucial recalculation step discussed in Sec. II. As the figure shows, without this step, the algorithm fails to find the division of the network into the two known groups. Furthermore, the modularity doesn’t reach nearly such high values as in the first two panels, indicating that the divisions suggested are much poorer than in the cases with the recalculation.

### V.3 Collaboration network

For our next example, we look at a collaboration network of scientists. Figure 10a shows the largest component of a network of collaborations between physicists who conduct research on networks. (The authors of the present paper, for instance, are among the nodes in this network.) This network (which appeared previously in Ref. PN03, ) was constructed by taking names of authors appearing in the lengthy bibliography of Ref. Newman03d,  and cross-referencing with the Physics E-print Archive at arxiv.org, specifically the condensed matter section of the archive where, for historical reasons, most papers on networks have appeared. Authors appearing in both were added to the network as vertices, and edges between them indicate coauthorship of one or more papers appearing in the archive. Thus the collaborative ties represented in the figure are not limited to papers on topics concerning networks—we were interested primarily in whether people know one another, and collaboration on any topic is a reasonable indicator of acquaintance.

The network as presented in Fig. 10a is difficult to interpret. Given the names of the scientists, a knowledgeable reader with too much time on their hands could, no doubt, pick out known groupings, for instance at particular institutions, from the general confusion. But were this a network about which we had no a priori knowledge, we would be hard pressed to understand its underlying structure.

Applying the shortest-path version of our algorithm to this network we find that the modularity, Eq. (5), has a strong peak at 13 communities with a value of $`Q=0.72\pm 0.02`$. Extracting the communities from the corresponding dendrogram, we have indicated them with colors in Fig. 10b. The knowledgeable reader will again be able to discern known groups of scientists in this rendering, and more easily now with the help of the colors. Still, however, the structure of the network as a whole and the of the interactions between groups is quite unclear.

In Fig. 10c we have reduced the network to *only* the groups. In this panel, we have drawn each group as a circle, with size varying roughly with the number of individuals in the group. The lines between groups indicate collaborations between group members, with the thickness of the lines varying in proportion to the number of pairs of scientists who have collaborated. Now the overall structure of the network becomes easy to see. The network is centered around the middle group shown in cyan (which consists of researchers primarily in southern Europe), with a knot of inter-community collaborations going on between the groups on the lower right of the picture (mostly Boston University physicists and their intellectual descendants). Other groups (including the authors’ own) are arranged in various attitudes further out.

One of the problems created by the sudden availability in recent years of large network data sets has been our lack of tools for visualizing their structure Newman03d . In the early days of network analysis, particularly in the social sciences, it was usually enough simply to draw a picture of a network to see what was going on. Networks in those days had ten or twenty nodes, not 140 as here, or several billion as in the world wide web. We believe that methods like the one presented here, of using community structure algorithms to make a meaningful “coarse graining” of a network, thereby reducing its level of complexity to one that can be interpreted readily by the human eye, will be invaluable in helping us to understand the large-scale structure of these new network data.

### V.4 Other examples

In this section, we briefly describe example applications of our methods to three further networks. The first is a non-human social network, a network of dolphins, the second a network of fictional characters, and the third not a social network at all, but a network of web pages and the links between them.

In Fig. 11 we show the social network of a community of 62 bottlenose dolphins living in Doubtful Sound, New Zealand. The network was compiled by Lusseau Lusseau03b  from seven years of field studies of the dolphins, with ties between dolphin pairs being established by observation of statistically significant frequent association. The network splits naturally into two large groups, represented by the circles and squares in the figure, and the larger of the two also splits into four smaller subgroups, represented by the different shades. The modularity is $`Q=0.38\pm 0.08`$ for the split into two groups, and peaks at $`0.52\pm 0.03`$ when the subgroup splitting is included also.

The split into two groups appears to correspond to a known division of the dolphin community Lusseau03a . Lusseau reports that for a period of about two years during observation of the dolphins they separated into two groups along the lines found by our analysis, apparently because of the disappearance of individuals on the boundary between the groups. When some of these individuals later reappeared, the two halves of the network joined together once more. As Lusseau points out, developments of this kind illustrate that the dolphin network is not merely a scientific curiosity but, like human social networks, is closely tied to the evolution of the community. The subgroupings within the larger half of the network also seem to correspond to real divisions among the animals: the largest subgroup consists almost of entirely of females and the others almost entirely of males, and it is conjectured that the split between the male groups is governed by matrilineage (D. Lusseau, personal communication).

Figure 12 shows the community structure of the network of interactions between major characters in Victor Hugo’s sprawling novel of crime and redemption in post-restoration France, Les Misérables. Using the list of character appearances by scene compiled by Knuth Knuth93 , the network was constructed in which the vertices represent characters and an edge between two vertices represents co-appearance of the corresponding characters in one or more scenes. The optimal community split of the resulting graph has a strong modularity of $`Q=0.54\pm 0.02`$, and gives 11 communities as shown in the figure. The communities clearly reflect the subplot structure of the book: unsurprisingly, the protagonist Jean Valjean and his nemesis, the police officer Javert, are central to the network and form the hubs of communities composed of their respective adherents. Other subplots centered on Marius, Cosette, Fantine, and the bishop Myriel are also picked out.

Finally, as an example of the application of our method to a non-social network, we have looked at a web graph—a network in which the vertices and edges represent web pages and the links between them. The graph in question represents 180 pages from the web site of a large corporation 555The graph is one of the test graphs from the graph drawing competition held in conjunction with the Symposium on Graph Drawing, Berkeley, California, September 18–20, 1996.. Figure 13 shows the network and the communities found in it by the shortest-path version of our algorithm. This network has one of the strongest modularity values of the examples studied here, at $`Q=0.65\pm 0.02`$. The links between web pages are directed, as indicated by the arrows in the figure, but, as discussed in Sec. III.1, for the purposes of finding the communities we ignore direction and treat the network as undirected.

Certainly it might be useful to know the communities in a web network; an algorithm that can pick out communities could reveal which pages cover related topics or the social structure of links between pages maintained by different individuals. Ideas along these lines have been pursued by, for example, Flake et al. FLGC02  and Adamic and Adar AA01 .

## VI Conclusions

In this paper we have described a new class of algorithms for performing network clustering, the task of extracting the natural community structure from networks of vertices and edges. This is a problem long studied in computer science, applied mathematics, and the social sciences, but it has lacked a satisfactory solution. We believe the methods described here give such a solution. They are simple, intuitive, and demonstrably give excellent results on networks for which we know the community structure ahead of time. Our methods are defined by two crucial features. First, we use a “divisive” technique which iteratively removes edges from the network, thereby breaking it up in communities. The edges to be removed are identified by using one of a set of edge betweenness measures, of which the simplest is a generalization to edges of the standard shortest-path betweenness of Freeman. Second, our algorithms include a recalculation step in which betweenness scores are re-evaluated after the removal of every edge. This step, which was missing from previous algorithms, turns out to be of primary importance to the success of ours. Without it, the algorithms fail miserably at even the simplest clustering tasks.

We have demonstrated the efficacy and utility of our methods with a number of examples. We have shown that our algorithms can reliably and sensitively extract community structure from artificially generated networks with known communities. We have also applied them to real-world networks with known community structure and again they extract that structure without difficulty. And we have given examples of how our algorithms can be used to analyze networks whose structure is otherwise difficult to comprehend. The networks studied include a collaboration network of scientists, in which our methods allow us to generate schematic depictions of the overall structure of the network and collaborations taking place within and between communities, other social networks of people and of animals, and a network of links between pages on a corporate web site.

The primary remaining difficulty with our algorithms is the relatively high computational demands they make. The fastest of them, the one based on shortest-path betweenness, operates in $`\mathrm{O}(n^3)`$ time on a sparse graph, which makes it usable for networks up to about $`\mathrm{10\hspace{0.17em}000}`$ vertices, but for larger systems it becomes intractable. Although the ever-improving speed of computers will certainly raise this limit in coming years, it would be more satisfactory if a faster version of the method could be discovered. One possibility is parallelization: the betweenness calculation involves a sum over source vertices and the elements of that sum can be distributed over different processors, making the calculation trivially parallelizable on a distributed-memory machine. However, a better approach would be to find some improvement in the algorithm itself to decrease its computational complexity.

Since the publication of our first paper on this topic GN02 , several other authors have made use of the shortest-path version of our algorithm. Holme et al. HHJ02  have applied it to a number of metabolic networks for different organisms, finding communities that correspond to functional units within the networks, while Wilkinson and Huberman WH03  have applied it to a network of relations between genes, as established by co-occurrence of names of genes in published research articles. An interesting application to social networks is the study by Gleiser and Danon GD04  of the collaboration network of early jazz musicians. They found, among other things, that the network split into two communities along lines of race, black musicians in one group, white musicians in the other. Guimerà et al. Guimera03  have applied the method to a network of email messages passing between users at a university, and found communities that reflect both formal and informal levels of organization. Tyler et al. TWH03  have also applied the algorithm to an email network, in their case at a large company, finding that the resulting communities correspond closely to organizational units. The latter work is interesting also in that it suggests a method for improving the speed of the algorithm: Tyler et al. calculate betweenness for only a subset, randomly chosen, of possible source vertices in the network, rather than summing over all sources. The size of the subset is decided on the fly, by sampling source vertices until the betweenness of at least one edge in the network exceeds a predetermined threshold. This technique reduces the running time of the calculation considerably, although the resulting estimate of betweenness necessarily suffers from the statistical fluctuations inherent in random sampling methods. This idea, or a variation of it, might provide a solution to the problems mentioned above of the high computational demands of our algorithms.

We are of course delighted to see our methods applied to such a variety of problems. Combined with the new algorithms and measures described in this paper, we hope to see many more applications in the future.

###### Acknowledgements.

The authors thank Steven Borgatti, Ulrik Brandes, Linton Freeman, David Lusseau, Mason Porter, and Douglas White for useful comments. Thanks also to Oliver Boisseau, Patti Haase, David Lusseau, and Karsten Schneider for providing the data for the dolphin network and to Douglas White for the karate club data. This work was funded in part by the National Science Foundation under grant number DMS–0234188 and by the Santa Fe Institute.