Relaxation Optimized Transfer of Spin Order in Ising Spin Chains

Relaxation optimized transfer of spin order in Ising spin chains

in
o
o

c

in

(N

>

in :
o .
in ■
o ■

^ :

Qk
-i— > ■
G ■
cd

3 :
cr

13

Dionisis Stefanatos^'Q Steffen J. Glaser, 2 and Navin Khaneja 1

'Division of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA
Department of Chemistry, Technische Universitat Munchen, 85747 Garching, Germany

(Dated: February 1, 2008)

In this manuscript, we present relaxation optimized methods for transfer of bilinear spin correla-
tions along Ising spin chains. These relaxation optimized methods can be used as a building block
for transfer of polarization between distant spins on a spin chain. Compared to standard techniques,
significant reduction in relaxation losses is achieved by these optimized methods when transverse
relaxation rates are much larger than the longitudinal relaxation rates and comparable to couplings
between spins. We derive an upper bound on the efficiency of transfer of spin order along a chain
of spins in the presence of relaxation and show that this bound can be approached by relaxation
optimized pulse sequences presented in the paper.

PACS numbers: 03.67.-a, 03.65.Yz, 82.56.Jn, 82.56.Fk

I. INTRODUCTION

Relaxation (dissipation and decoherence) is a charac-
teristic feature of open quantum systems. In practice,
relaxation results in loss of signal and information and ul-
timately limits the range of applications. Recent work in
optimal control of spin dynamics in the presence of relax-
ation has shown that these losses can be significantly re-
duced by exploiting the structure of relaxation 0,000].
This has resulted in significant improvement in sensitiv-
ity of many well established experiments in high resolu-
tion nuclear magnetic resonance (NMR) spectroscopy. In
particular, by use of optimal control methods, analytical
bounds have been achieved on the maximum polarization
or coherence that can be transferred between coupled
spins in the presence of very general decoherence mecha-
nisms. In this paper, we look at the more general problem
of transfer of coherence or polarization between distant
spins on an Ising spin chain in the presence of relax-
ation. This problem is ubiquitous in multi-dimensional
NMR spectroscopy, where polarization is transferred be-
tween distant spins on a chain of coupled spins. Spin
(or pseudo spin) chains also appear in many proposed
quantum information processing architectures |3, |6| .

The system that we study in this paper is a linear chain
of n weakly interacting spins 1/2 placed in a static exter-
nal magnetic field in the z direction (NMR experimental
setup), with Ising type couplings of equal strength be-
tween nearest neighbors, see Fig. ^ The free evolution
Hamiltonian of the system has the form

27rJ hzl{i+i)z,

i=l

where u>i is the Larmor frequency of spin i and J is the
strength of the coupling between the spins. In a suitably
chosen (multiple) rotating frame, which rotates with each

OrO

J

J

J

J

n-1

FIG. 1: The system that we study in this paper is a linear
chain of n weakly interacting spins 1/2 with Ising coupling
between next neighbors. The coupling constant is the same
for all pairs of connected spins.

spin at its resonance (Larmor) frequency, the free evolu-
tion Hamiltonian simplifies to

n-i
i=i

(1)

Motivated by NMR spectroscopy of large molecules
in solution, we assume that the relaxation rates of the
longitudinal operators with components only in the z
direction, like Ii z , 2Ii Z Iu+i-\ z is negligible compared
with relaxation rates for transverse operators like Ii X ,
2IizI(i+i)y 0- The transverse relaxation is modeled by
the Lindbladian with the general form

L(P)

E

it a,

[Iiz[Iiz,p]] +X) 7f&

ij {21izljz[21izljz, p]] ■

'Electronic address: stefanat@fas.harvard.edu

In liquid state NMR spectroscopy, the two terms of
the Lindblad operator model the relaxation mechanism
caused by chemical shift anisotropy and dipole-dipole in-
teraction respectively 0. Here we neglect any interfer-
ence effects between these two relaxation mechanisms 0] .
Relaxation rates , bij depend on various physical pa-
rameters, such as the gyromagnetic ratios of the spins,
the internuclear distance, the correlation time of molecu-
lar tumbling etc. We define the net transverse relaxation
rate for spin i as k l a — fli+J^j Without loss of gen-
erality in the subsequent analysis, we assume that k l a are
equal and we denote this common transverse relaxation
rate by k.

The time evolution (in the rotating frame) of the spin
system density matrix p is given by the master equation

P = -i[H,p] + L(p)

(2)

where H = H c + H r f and H r f is the control Hamiltonian.
In the NMR context, the available controls are the com-
ponents of the transverse radio-frequency (RF) magnetic
field. It is assumed that the resonance frequencies of the
spins are well separated, so that each spin can be selec-
tively excited (addressed) by an appropriate choice of the
components of the RF field at its resonance frequency.

Consider now the problem of optimizing the polariza-
tion transfer

Iu

(3)

along the linear spin chain shown in Fig. ^ in the
presence of the relaxation mechanisms mentioned above.
This problem can be stated as follows: Find the opti-
mal transverse RF magnetic field in the control Hamil-
tonian H r f such that starting from p(0) = I\ z and
evolving under Eq. J5J) the target expectation value
(Inz) — kr{p(t)I nz } is maximized.
To fix ideas, we analyze the case when n = 3

Iu

hz

(4)

This transfer is achieved conventionally using INEPT like
pulse sequences 0,^,^3- Under the conventional trans-
fer method, the initial state of the system I\ z evolves
through the following stages

hz — * 2-Zi z /2

2I 2z hz

'3z-

(5)

In the first stage of the transfer, I\ z — > 2Ii z I 2zi spin 3
is decoupled from the chain using standard decoupling
methods [12| and the initial polarization I\ z on spin 1 is
rotated by RF field (an appropriate 7r/2 pulse) to coher-
ence l\ x , which then evolves under coupling Hamiltonian
2I\ z Iiz to 21\ y I 2z . When the expectation value (2I\ y Ii z )
is maximized, another 7r/2 pulse is used to rotate 2I\ v l 2z
to 2/i 2 /22- This is the INEPT pulse sequence. The next
stage, 2I lz I 2z — * 2I 2z I 3z , is the so-called spin order trans-
fer. Fig. |3 shows the population inversion corresponding
to this transfer. By a suitable ir/2 rotation of spin 2,
the density operator 2I lz I 2z is transformed to 2I lz I 2x ,
which then evolves to 2I 2x I 3z . When the expectation
value (2I 2x I 3z ) is maximized, another ir/2 pulse is used
to rotate 2I 2x I 3z to the spin order 2I 2z I 3z . This is the
Concatenated INEPT (CINEPT) pulse sequence. The fi-
nal stage transfer, 2I 2z I 3z — > I 3z , is similar to that in the
first stage and is accomplished by the INEPT pulse se-
quence. The efficiency of these transfers is limited by the
decay of transverse operators, due to the phenomenon of
relaxation.

In our recent work on relaxation optimized control of
coupled spin dynamics , we showed that the efficiency
of the first and the last step (the two INEPT stages)
in Eq. (jSJ) can be significantly improved by controlling

■PPa

■ aaa

FIG. 2: Energy level diagram for the spin order transfer
2Ii z l2z — » 2l2 Z l3z- The dark circles represent population ex-
cess. State a corresponds to spin up and state /3 to spin down.

precisely the way in which magnetization is transferred
from longitudinal operators to transverse operators. In
other words, instead of using 7r/2 (hard) pulses to rotate
longitudinal to transverse operators (and the inverse), we
can exploit the fact that longitudinal operators are long
lived by making these rotations gradually, saving this way
magnetization. This transfer strategy is called relaxation
optimized pulse element (ROPE).

In this article we derive relaxation optimized pulse
sequences for the intermediate transfer (the CINEPT
stage),

2/i z /2

2hJ,

(6)

This relaxation optimized transfer of spin order can then
be used as a building block for the polarization transfer
© through the scheme

Iiz ~ * 2I\ Z I 2

2I,J:

2z*3z

21,

In

(7)

II. THE OPTIMAL CONTROL PROBLEM AND
AN UPPER BOUND FOR THE EFFICIENCY

In this section, we formulate the problem of trans-
fer in Eq. © as a problem of optimal control and
derive an upper bound on the transfer efficiency. To
simplify notation, we introduce the following symbols
for the expectation values of operators that play a part
in the transfer. Let z\ = (2I\ z I 2z ), x\ = (2Ii z I 2x ),
y 2 = (V2{2I lz I 2y I 3z + I 2y /2)), x 3 = -(2I 2x I 3z ) and
23 = (2I 2z I 3z ). As a control variable we use the trans-
verse RF magnetic field, pointing say in the y direction
(in the rotating frame), so H r f = ui y (t)I 2y . Note that
u! y (t) is the component of the field in the rotating frame,
so it is actually the envelope of the RF field. The carrier
frequency of the RF field is the resonance frequency of
spin 2. Using Eq. @, we find that the evolution of the

3

system, in time units of 1/ , is given by

Zl '

"

Qy

±1

-£

-1

m

1

-€

-1

±3

1

-€

f2 a

" Zl '

X\

2/2

^3

. Z3 .

(8)

where fi„(i) = w v (t)/(irJV2) and £ = k/{jV2). The
initial condition is (zi, Xi,y2, x 3 , 23) — (1, 0, 0, 0, 0).

The efficiency of the conventional method (CINEPT)
for transfer Z\ — > Z3, can be easily found. At t — 0, z\ is
transferred to X\ by application of a (tt/2)^ pulse on spin
2. Couplings evolve xi to j/2 which further evolves to X3.
As a function of time, 2:3 (t) = e - ^' sin 2 (t/v / 2)- This is
maximized for t m = v / 2cot~ 1 (^/v / 2)- At i = i m a second
{tt/2) v pulse is applied on spin 2, rotating the maximum
value from 23 to 23. This value is the efficiency 77c"/ of
the conventional method

Vci

= exp(-^V2cot- 1 (^/v / 2))sin 2 (cot- 1 (,e/v / 2)) . (9)

A better efficiency can be achieved if we store magne-
tization in the decoherence free longitudinal operators Zi
while the system is evolving. This is done by rotating z%
to xi gradually, instead of using tt/2 hard pulses. This
is the physical concept behind the relaxation optimized
transfer strategy. For the specific transfer examined in
this article, we first find an upper bound for the maxi-
mum efficiency and in the next section we calculate nu-
merically the magnetic field £l y (t) that approaches this
bound.

In order to derive the upper bound, we use an aug-
mented system instead of the original one iJHJ. The aug-
mentation is done in two steps. First, we suppose that
we can rotate z\ to X\ and X3 to Z3 independently using
two different controls, say Qi(t) and SI3 (i) , instead of the
common control Q y (t). Next, we provide 7/2 with a relax-
ation free partner zi and with a control Q2 (i) which can
rotate yi to z%. The augmented system is

" Zl '

"

-Oi

-

" Zl '

Zl

-fia

Z2

Z3

Z3

±1

ill

-e

-1

Xl

2/2

a 2

t

-£

-f

2/2

. £3 .

.

<h

1

-€ -

. X3 .

(10)

Observe that system (fTU)) reduces to system © for fii =
— O3 = fly and CI2 — 0. Thus, if we know the maximum
achievable value of z 3 starting from (1,0,0,0,0,0) and
evolving under system l|10|l , then this is an upper bound
for the maximum achievable value of Z3 with evolution
described by the original system JSJ.

and

Let n

=

+ zf,

T2

= vvi + z

'■3

\J x\ + z\. Using fli(t) we can control the angles 9i,
shown in Fig. [21 independently. If we assume that the
control can be done arbitrarily fast as compared to the
evolution of couplings or relaxation rates then we can

z il r z a

z A r

1 J 2

FIG. 3: Auxiliary variables rt.

think of 9i as control variables. The equations for the
evolution of r, are

fl

r2

U1U2

h 2
U2U3

-u 2 u 3

ri

T2

. r3 .

(11)

where the new control parameters are Ui(t) = cos(9i(t)).
The goal is to find the largest achievable value of T3 start-
ing from (ri,r2,r3) = (1,0,0) by appropriate choice of
Ui{t). This problem can be solved analytically. The opti-
mal solution is characterized by maintaining vanishingly
small values of dr2/dt, i.e., (uiri — u$rs)/u2T2 = £ and
by U3T3/uiri = k, where

(12)

The maximum achievable value of T3 is also k. We prove
it in the following.

Using variables p t = rf, m, = uin/ \Jy^ (u.n) 2 , and

dr/dt = (uiTi) 2 , equation (jll|) can be re-written as

d_

dr

where

Pi( T )
M T )

A :

diag{Arn{r)m r {t))

-i -1
1 -i -1

1 -£

(13)

(14)

m T = (mi,TO2,m3) and diag(X) represents the vector
containing diagonal entries of the square matrix X. The
goal is to find the controls TOj(r) m f( T ) — 1) and the
largest achievable value of P3 starting from {pi,P2,P3) —
(1,0,0). Eq. O implies that

diag(A / m(r)m T \r)dr)
Jo

Let M — m(T)m T (r)dr . Note that M is a symmet-
ric, positive semidefinite matrix. By definition Pi(r) >
for all t. At final time T, we must have Pi(T) = and
P2{T) = 0, as any nonzero value of Pi(T) or piiT) can

'MT)

Mo)

Mo)

Mo)

4

be partly transferred to and the final value of Ps{T)
further increased. Since (pi(0),p 2 (0),P3(0)) = (1,0,0),
it implies that M should be such that [AM)n = — 1,
(AM)%2 — and (AM)^ is maximized over all positive
scmidefinite M satisfying the above constraints. This
problem is a special case of a semidcfmitc programming
problem |13| . If the symmetric part of matrix A is neg-
ative definite, as in our case l|14|) . then it can be shown
that the optimal solution M to the above stated semi-
definite programming problem is a rank one matrix |l4j ,
i.e., M = mm T for some constant m and therefore the
ratio Uiri/ujTj in (|llfl is constant throughout. The con-
dition (AM)22 = 0, implies (u\ri — u^r^)/^ = U2T2- Sub-
stituting for U2T2 in JUJ, we obtain

h

T\

h _

u\uiii

r3

(15)

We now simply need to maximize the gain in r% to loss
in r%, i.e., the ratio r^/i—ri). This yields u^r^/uiri = k,
with k given in O- The corresponding maximum effi-
ciency for transfer rj — > is also ft. This is the maximum
efficiency for transfer z\ — > Z3 under the augmented sys-
tem l|10|). and thus an upper bound for the efficiency of
the same transfer under the original system

TABLE I: For various values of £ £ [0, 1], the optimal values of
A, a and the corresponding efficiency are shown. We present
also for comparison the efficiency achieved by the steepest
descent method.

I

A

a

Gaussian Pulse

Steepest Descent

1.00

1.11

1.30

0.2510

0.2512

0.95

1.09

1.32

0.2661

0.2662

0.90

1.07

1.34

0.2824

0.2825

0.85

1.05

1.36

0.3000

0.3001

0.80

1.03

1.38

0.3190

0.3191

0.75

1.02

1.39

0.3396

0.3397

0.70

1.00

1.41

0.3619

0.3620

0.65

0.98

1.43

0.3861

0.3863

0.60

0.97

1.44

0.4124

0.4126

0.55

0.96

1.44

0.4410

0.4413

0.50

0.95

1.44

0.4721

0.4726

0.45

0.94

1.45

0.5060

0.5067

0.40

0.93

1.46

0.5428

0.5439

0.35

0.92

1.46

0.5830

0.5846

0.30

0.91

1.46

0.6270

0.6292

0.25

0.90

1.47

0.6750

0.6780

0.20

0.89

1.48

0.7277

0.7315

0.15

0.88

1.48

0.7855

0.7900

0.10

0.85

1.52

0.8494

0.8536

0.05

0.79

1.60

0.9203

0.9232

0.00

0.73

1.71

0.9999

1.0000

III. NUMERICAL CALCULATION OF THE
OPTIMAL RF FIELD AND DISCUSSION

Having established an analytical upper bound (|12fl for
the efficiency, we now try to find numerically a RF field
£l y (t) that approaches this bound, for each value of the
parameter £. We emphasize that in this section the orig-
inal system JSJ is employed.

At first, we use a numerical optimization method based
on a steepest descent algorithm. For the application of
the method we use a finite time window T. For values of
normalized relaxation £ in (0—1), a time interval T =
10 (normalized time units) is enough. For larger values
of £ we can use even shorter T. The optimal RF field
Q y {t) that we find with this method, for various £, is
shown in Fig. 0] Note that as £ increases, the optimal
pulse becomes shorter in time and acquires a larger peak
value. The reason for this is that for larger £ the transfer
z\ — > Z3 should be done faster, in order to reduce the time
spent in the transverse plane and hence the relaxation
losses. Now observe that the optimal pulse shape can be
very well approximated by a Gaussian profile of the form

fl y (t) — Aexp

t-T/2
V2a

(16)

with A, a appropriately chosen. As a result, the effi-
ciency that we find using the appropriate Gaussian pulse
is very close to that we find using the original pulse. This
suggests that instead of using the initial numerical opti-
mization method, we can use Gaussian pulses of the form

Hlfty . optimized with respect to A and a for each value of
£. The optimal A, a are found by numerical simulations.
For each £ we simulate the equations of system (jHJ) with
Fly(t) given by Eq. (|16fl . for many values of A and a.
We choose those values that give the maximum 23 (T).
In table I, we show the optimal A, a for various values
£ £ [0, 1]. We also show the corresponding efficiency, as
well as the efficiency achieved by the initial numerical
optimization method. Observe how close lie these two
groups of values. The choice of the Gaussian shape is
indeed successful.

Fig. 13 shows the efficiency of the conventional method
(CINEPT), i.e., i] C i from Eq. JjJ), the efficiency of our
method (SPORTS ROPE, SPin ORder TranSfer with Re-
laxation Optimized Pulse Element) , and the upper bound
k from Eq. I|12l) , for the values of relaxation parameter £
shown in table I. Note that for large £ (large relaxation
rates), SPORTS ROPE gives a significant improvement
over CINEPT. Also note that it approaches fairly well
the upper bound.

Using the Gaussian pulse shape we can get a quantita-
tive impression of the robustness of SPORTS ROPE. In
Fig. [S] we give a gray-scale topographic plot of the effi-
ciency (2:3 (T) with T = 10) as a function of A and a for
£ = 1. The maximum value can be found from table I and
is 0.2510. The white region corresponds to values > 0.24,
while the black region to values < r\ci{£, = 1) = 0.1727.
The intermediate gray regions correspond to values be-
tween these two limits. Obviously, SPORTS ROPE is
quite robust.

In Fig. Ufa), we plot the time evolution of the various

5

02468 10 02468

(a) t (b) t

02468 10 02468 10

(c) t (d) t

1 2 3 4 5 0.5 1 1.5 2 2.5

(e) t (f) t

FIG. 4: Optimal pulse (dashed line) calculated using a numerical optimization method based on a steepest descent algorithm,
for various values of the normalized relaxation parameter £. The Gaussian pulse (solid line) approximates very well the optimal
pulse shape and gives a similar efficiency. This suggests that instead of using the initial numerical optimization method, we
can use Gaussian pulses of the form I16H . optimized with respect to A and a for each value of £.

transfer functions (expectation values of operators) that
participate in the transfer z\ — > z 3 , when the optimal
Gaussian pulse for £ = 1, shown in Fig. Hfc), i s applied to
system ©. Observe the gradual building of the interme-
diate variables x±, y% and X3. Note that dr^jdt — 1/2 7^ 0.
There is no contradiction with the optimality condition
dr 2 /dt = derived in section [n] during the calculation of

the upper bound, since this condition refers to the aug-
mented system (I1U|I and not the original one JSJ used
here. In Fig. Efb) we plot the angle #3 = tan -1 (Z3/X3)
of the vector r% with the x axis, as a function of time.
Observe that initially r 3 is parallel to x axis (8 3 = 0), but
under the action of the Gaussian pulse is rotated gradu-
ally to z axis (#3 = 7r/2). This gradual rotation of r 3 (as

6

0.8

>>

0.6

o

'o

0.4

0.2

Q

» 1 r

g
o

5

° S

□

x □

o CINEPT

« SPORTS ROPE

" UPPER BOUND

o * □
° * □

»

o

□

« □

x □
x □
O x □

° „ ' « □ n
O x □ _

0.2

0.4 0.6

0.8

FIG. 5: Efficiency for the conventional method (CINEPT),
Eq. and for our method (SPORTS ROPE), for the values
of £ shown in table I. The upper bound (1121 for the efficiency
is also shown.

well as of n) is a characteristic feature of the SPORTS
ROPE transfer scheme.

We remark that for the general transfer I\ z — >
I nz , more than one intermediate steps 2Iu^x)zhz ~ *
2IizI(i+i) z are necessary. Since the equations that de-
scribe the i th transfer are the same as (jHJ, we just need
to apply the same Gaussian pulse but centered, in the
frequency domain, at the resonance frequency of spin i.
In this sequence of Gaussian pulses we should add at the
beginning and at the end the optimal pulses for the first
and the final step, respectively, see Eq. (J7J). These pulses
can be calculated using the theory presented in Fi-
nally, note that the same scheme can be used for the
coherence transfer I± a — > I n p, where a,/3 can be x or y.
We just need to add the initial and final tt/2 pulses that

accomplish the rotations I\ a — > I\ z , I n

IV. CONCLUSION

0.4 0.8 1.2 1.6

In this paper, we derived an upper bound on the ef-
ficiency of spin order transfer along an Ising spin chain,
in the presence of relaxation, and calculated numerically
relaxation optimized pulse sequences approaching this
bound. Using these methods, a significant reduction in
relaxation losses is achieved, compared to standard tech-
niques, when transverse relaxation rates are much larger
than the longitudinal relaxation rates and comparable
to couplings between spins. These relaxation optimized
methods can be used as a building block for transfer of
polarization or coherence between distant spins on a spin
chain. This problem is ubiquitous in multi-dimensional
NMR spectroscopy and is also interesting in the context
of quantum information processing.

FIG. 6: Gray-scale topographic plot of the efficiency
as a function of the amplitude A and the standard
deviation a of the Gaussian pulse for f = 1. The max-
imum value can be found from table I and is 0.2510.
The white region corresponds to values > 0.24, while
the black region to values < r)ci(£=l) = 0.1727. The
intermediate gray regions correspond to the intervals
[0.23, 0.24), [0.22, 0.23), [0.21, 0.22), [0.20, 0.21), [0.1727, 0.20)
(from white to black).

Acknowledgments

N.K. acknowledges Grants AFOSR FA9550-04-1-0427,
NSF 0133673 and NSF 0218411. S.J.G. thanks the
Deutsche Forschungsgemeinschaft for Grant Gl 203/4-2.

[1] N. Khaneja, T. Reiss, B. Luy, and S.J. Glaser, J. Magn.

Reson. 162, 311 (2003).
[2] N. Khaneja, B. Luy, and S.J. Glaser (2003) Proc. Natl.

Acad. Sci. U.S.A 100, 13162 (2003).
[3] D. Stefanatos, N. Khaneja, and S.J. Glaser, Phys. Rev.

A, 69, 022319 (2004).
[4] N. Khaneja, Jr. Shin Li, C. Kehlet, B. Luy, S.J. Glaser,

Proc. Natl. Acad. Sci. USA. 101, 14742-47 (2004).
[5] B.E. Kane, Nature 393, 133 (1998).
[6] F. Yamaguchi, Y. Yamamoto, Appl. Phys. A 68 (1999).
[7] J. Cavanagh, W.J. Fairbrother, A.G. Palmer III, and N.J.

Skelton, Protein NMR Spectroscopy (Academic Press,

New York, 1996).
[8] M. Goldman, Interference effects in the relaxation of a

7

1

2 0.8
_o

|0.6
=H

£0.4
a

S 0.2

z l

x 1 \

(a)

t

10

(b)

FIG. 7: (a) Time evolution of the various transfer functions (expectation values of operators) participating in the transfer zi — >
Z3, when the optimal Gaussian pulse for £ = 1, shown in Fig. |IJc), is applied to system JSJ. (b) The angle 83 — tan -1 (23/0:3)
of the vector 7-3 with the X cLXlS, ELS £L function of time. Observe that initially r-j, is parallel to x axis (63 = 0), but under the
action of the Gaussian pulse is rotated gradually to z axis (63 — tt/2).

pair of unlike spin-1/2 nuclei, J. Magn. Reson. 60, 437
(1984).

[9] G. A. Morris, R. Freeman, J. Am. Chem. Soc. 101, 760

(1979) .

[10] D. P. Burum, R. R. Ernst, J. Magn. Reson. 39, 163

(1980) .

[11] A. Majumdar, E. P. Zuiderweg, J. Magn. Reson. A 113,
19-31 (1995).

[12] R.R. Ernst, G. Bodenhausen, A. Wokaun, Principles of
Nuclear Magnetic Resonance in One and Two Dimen-
sions, Clarendon Press, Oxford, 1987.

[13] L. Vandenberghe, S. Boyd, SIAM Review 38, 49-95
(1996).

[14] D. Stefanatos and N. Khaneja, \math. OC /0504308
(2005).