text:
# Mid-Infrared Imaging and Modelling of the Dust Shell around Post-AGB star HD 187885 (IRAS 19500-1709)

## 1. Introduction

The circumstellar envelopes (CSEs) of post-AGB stars are cool and dust rich and hence radiate most of their energy at mid-infrared wavelengths. IRAS 19500-1709 is associated with the high galactic latitude F2-3I (Parthasarathy, Pottasch & Wamsteker 1988) post-AGB star HD 187885 and has the double-peaked spectral energy distribution (SED) typical of post-AGB stars with detached CSEs (Hrivnak, Kwok & Volk 1989). The expansion velocity of the envelope, based on CO-line emission, is 11 km s-1 with wings up to 30 km s-1 (Likkel et al. 1987). The envelope is carbon rich, showing no OH or H2O maser emission (Likkel 1989). It has a broad emission feature from 10-13$`\mu `$m peaking at 12$`\mu `$m which could be attributed to a polycyclic aromatic hydrocarbon (PAH), such as chrysene (Justtanont et al. 1996). It has a weak 21$`\mu `$m feature (Justtanont et al. 1996) and a broad feature around 30$`\mu `$m which has been modelled with non-spherical magnesium sulphide (MgS) dust grains (Hony, Waters & Tielens 2002). Imaging polarimetry at near-infrared wavelengths shows that IRAS 19500-1709 has a bipolar structure in scattered light (Gledhill et al. 2001). We present mid-infrared imaging of IRAS 19500-1709 using the OSCIR camera mounted on the 8.1-m Gemini North Telescope which provides information on the inner part of the CSE. We fit the spectral energy distribution (SED) from UV to sub-mm wavelengths and simulate the images using a 2-D dust radiation transport (RT) code in order to derive the physical and chemical properties of the dust in the CSE.

## 2. Observations

Using the OSCIR camera on Gemini North, we obtained 10 and 20$`\mu `$m (N and Q3 band) images of IRAS 19500-1709 with spatial resolution of 0.5 arcsecs (N band) to gain information on the inner part of the CSE (Fig.1). The object is extended relative to the point spread function (PSF).

The N band suggests the object is elongated in a N-S direction. The outer contours appear rounder in the Q band image but there is evidence of elongation in the centre in a N-S direction.. Also the brightness peak in the N-band image appears approximately 0.3 arcsecs off centre to the east, relative to the outer contours.

## 3. Modelling of the Dust Shell

We modelled the images and SED of IRAS 19500-1709 using a 2-D axisymmetric RT dust code (Efstathiou & Rowan-Robinson 1990) to derive the physical and chemical properties of dust in the envelope. The SED is constrained by the available observational data, for IRAS 19500-1709, taken from the literature which includes a 7.6-23.6$`\mu `$m UKIRT spectrum (CGS3) provided by K.Justannont, optical and near-IR photometry (Hrivnak, Kwok & Volk 1989), 8.5, 10 and 12.2$`\mu `$m flux (Meixner et al. 1997), IRAS fluxes, IUE data111Based on INES data from the IUE satellite, ISO SWS222ISO SWS06 observing program rszczerba-PPN30 and LWS333ISO LWS01 observing program mbarlow-dust 3 spectra and SCUBA flux (Van der Veen et al. 1994). Fig. 2 shows our best fit to date.

In this model we use a spherically symmetric shell with dust composed of 0.01-2$`\mu `$m amorphous carbon (amC), silicon carbide (SiC) and MgS grains. The amC dust fits the continuum of the SED quite well. We include SiC to account for the presence of the infrared emission feature from 10-13$`\mu `$m. We attempted to fit this feature using PAHs but found that they produced emission features in the 16-25$`\mu `$m region of the spectrum which are not seen in the data. The MgS is included as being the possible carrier of the 30$`\mu `$m band, although it is possible that this feature is due to an unknown carbonaceous component. We used optical constants for MgxFe(1-x)S from Begemann et al ((1994) with x=0.9 as an approximation to MgS (also see Hony et al. 2003). These data do not extend below 10$`\mu `$m however, so it was necessary to approximate the optical and UV extinction cross-sections to ensure that the MgS grains absorb sufficient radiation at short wavelengths to radiate in the mid-IR. This was done by adopting the optical constants of graphite below 4$`\mu `$m and interpolating between 4 and 10$`\mu `$m to produce a smooth extinction curve. We assume all dust grains are spherical with a size distribution of n(a)$`\alpha `$a-q with q=5. The dust is composed of 30% amC (temp. 87-167K), 40% SiC (temp. 79-150K) and 30% MgS (temp. 59-121K). The central star is assumed to be a black body with T=7500K. The ratio of the stellar radius to the inner radius of the shell is determined to be 9x10-5 in our model and Av=1.5. The ratio of the inner shell radius to the outer shell radius is 0.2 which produces sufficient far-IR emission to fit the data longward of 100$`\mu `$m. Extending the outer boundary any further than this produces too much far-IR emission.

In Fig.3a we show a raw model image corresponding to the best fit envelope parameters. In Fig.3b the model image is smoothed using a gaussian filter to simulate the OSCIR N-band image. The degree of smoothing (FWHM of the Gaussian) is determined by matching profiles between the model and OSCIR images. A gaussian smoothing function with a FWHM of 40 pixels provides the best match. This corresponds to a model pixel size of 0.01 arcsec. Using this procedure we estimate that the inner radius of the detached shell is 0.2 arcsec.

## 4. Future Work

There appears to be excess emission at 7-10$`\mu `$m as well as the broad feature around 30$`\mu `$m. We have attempted to fit this but the MgS peak is too narrow. This may be due to the MgS grains being too cold or because we have used spherical grains. Variations in grain shape can influence the emission profile with elliptical grains causing a wider peak (Hony, Waters & Tielens 2002). We will be running axisymmetric models to determine whether we can reproduce the elongation seen in the OSCIR images.

## References

Begemann, B. et al. 1994, ApJ, 423, L71

Efstathiou, A; & Rowan-Robinson, M. 1990, MNRAS, 245, 275

Gledhill, T. M. et al. 2001, MNRAS, 322, 321

Hony, S; Waters, L. B. F. M; & Tielens, A. G. G. M. 2002, A&A, 390, 533

Hony, S. et al. 2003, A&A, 402, 211

Hrivnak, B. J; Kwok, S; & Volk, K. M. 1989, ApJ, 346, 265

Justtanont, K. et al. 1996, A&A, 309, 612

Likkel, L. et al. 1987, A&A, 173, L11

Likkel, L. 1989, ApJ, 344, 350

Meixner, M. et al. 1997, ApJ, 482, 897

Parthasarathy, M; Pottasch, S. R; & Wamsteker, W. 1988, A&A, 203, 117

Van der Veen, W. E. C. J. et al. 1994, A&A, 285, 551


format: text

text:
# An Introduction to Using Software Tools for Automatic Differentiation (Revision A) ANL/MCS-TM-254. This work was supported by the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing Research, U.S. Department of Energy, under Contract W-31-109-ENG-38.

## 1 Intention and References

This document explains how to use the following software tools for automatic differentiation (AD) to generate first-order derivative code for a small example problem.

* ADIFOR 2.0 Revision D (http://www.mcs.anl.gov/adifor)
* ADIC 1.1 (http://www.mcs.anl.gov/adic)
* ADOL-C 1.8 (http://www.math.tu-dresden.de/wir/project/adolc)
* TAPENADE 1.0-alpha (http://www-sop.inria.fr/tropics)

The document does not describe the AD algorithms used by these software tools. Various publications cover the theoretical foundations of the algorithms referenced in our discussion. For an introduction and more advanced issues we point the reader to . Three workshops have been dedicated to AD, and the proceedings  contain numerous references to other AD-related papers that have appeared in scientific journals and proceedings of international conferences. A large variety of successful applications of AD software to real-world applications are discussed there as well.

## 2 Getting Started

For information on how to obtain, install, and use the AD tools discussed in this document, visit the Web sites listed in Section 1. The examples in this document are discussed in the context of the computer environment found on a laptop computer, specifically an Intel Pentium system running Mandrake Linux Version 8.1. ADIFOR, ADIC, and TAPENADE are installed in the directory /home/uwe/ADTOOLS. The ADOL-C library and include files are assumed to be in the same directory as the source files.

ADIFOR, ADIC, and TAPENADE are available on the Web and can be accessed online by using the corresponding Web servers. Visit

* http://www.mcs.anl.gov/autodiff/adiforserver for ADIFOR,
* http://www.mcs.anl.gov/autodiff/adicserver for ADIC, and
* http://tapenade.inria.fr:8080/tapenade/index.jsp for TAPENADE.

All source files presented in this document can be downloaded from

http://www-unix.mcs.anl.gov/̃naumann/ad\_tools.html.

The makefiles have to be adapted to the user’s environment.

## 3 Explosion Equation

We have selected a variant of the Bratu problem  to illustrate the use of source transformation AD tools. The code in Appendix A models the thermal explosion of solid fuels, which can be described by the system of differential equations

$$x^{\prime \prime }(\tau )+se^{\frac{x(\tau )}{1+tx(\tau )}}=0,$$

where $`\tau (1,1)`$ and $`x(1)=x(1)=0.`$ The problem has been discretized by using step size $`h`$ as

$$F_i=x_{i1}2x_i+x_{i+1}+h^2[f_{i1}+10f_i+f_{i+1}]/12$$

for $`i=1,\mathrm{},10000,`$ with $`x_0=x_{10001}=0`$ and $`f_i=s\mathrm{exp}(x_i/(1+tx_i)).`$ Of interest are the derivatives of the component functions $`F_i`$ with respect to the current state $`x_i`$ as well as the parameters $`s`$ and $`t.`$

We consider the following problems:

1. ADIFOR

1. Use the forward vector mode to compute the Jacobian of $`F`$ with respect to $`x`$ and the parameters $`s`$ and $`t`$. In our implementation, $`F`$ is represented by the variable f, $`x`$ is represented by x, and $`s`$ and $`t`$ are combined into the parameter vector prm. All three program variables are declared as arrays of double-precision floating-point numbers, namely,

    * double precision x(7), prm(2), f(dim) in Fortran and
        * double x, prm, f in C.

    (See Section 3.1.)
    2. Use Curtis-Powell-Reid  seeding to compute the compressed Jacobian in forward vector mode (see Section 3.1.1).
    3. Use the SparsLinC library to compute the sparse Jacobian by sparse forward mode (see Section 3.1.2).
2. ADOL-C

1. Use one of the easy-to-use drivers to compute the Jacobian of $`F`$ with respect to $`x`$ and the parameters $`s`$ and $`t`$.
    2. Compute the sparsity structure of the Jacobian.
    3. Use the low-level routines for forward and reverse mode to compute the Jacobian of $`F`$ with respect to $`x`$ and the parameters $`s`$ and $`t`$.

(See Section 3.3.)
3. ADIC

1. Compute the Jacobian of $`F`$ with respect to $`x`$ in forward vector mode (see Section 3.2).
4. TAPENADE

1. Compute the Jacobian of $`F`$ with respect to both $`x`$ and the parameters $`s`$ and $`t`$ in reverse mode (see Section 3.4).

### 3.1 ADIFOR 2.0

ADIFOR generates a differentiated version of the subroutine expl(...) with the following header.

```
      subroutine g_expl(g_p_, dim, parmax, x, g_x, ldg_x, prm, g_prm,
     + ldg_prm, f, g_f, ldg_f)
        integer dim, parmax
        double precision x(dim), prm(parmax)
        double precision f(dim)

        integer g_pmax_
        parameter (g_pmax_ = 9)
        integer g_p_, ldg_x, ldg_prm, ldg_f
        double precision g_x(ldg_x, dim), g_prm(ldg_prm, parmax),
     +    g_f(ldg_f, dim)

```

The parameter g\_pmax\_ is the maximal number of directional derivatives that can be computed. Its value is set by using AD\_PMAX in explosion.adf. The parameter is required because of the lack of dynamic memory allocation in Fortran 77. The derivative components of all scalar program variables are allocated as vectors of length g\_pmax\_. For example, in order to compute the Jacobian matrix of f with respect to x, the value of g\_pmax\_ must be at least equal to 7, so that the derivative components of x can accommodate the identity in $`IR^7.`$

The actual number of directional derivatives g\_p\_ must be less than or equal to g\_pmax\_. The parameter g\_p\_ determines the upper bound for the loops that compute the values of the derivative components.

To compute the Jacobian of f with respect to both x and prm, we set g\_pmax\_ to be greater than or equal to the sum of the numbers of elements in both vectors, that is, $`9=7+2=dim+\mathrm{𝚙𝚊𝚛𝚖𝚊𝚡}`$. Consequently, g\_x contains the first dim columns of the seed matrix and g\_prm its last parmax columns. The argument g\_f is the transpose of the $`(\mathrm{𝚍𝚒𝚖}\times \mathrm{𝚕𝚍𝚐}\mathrm{\_}𝚏)`$-Jacobian, where $`\mathrm{𝚕𝚍𝚐}\mathrm{\_}𝚏`$ must be initialized to 9. All this is done in the driver routine, which is shown in Appendix B.3.2.

The “dense:” section of the makefile in Appendix B.2 builds an executable that produces the following result.

```
-1.88  1.01  0.    0.    0.    0.    0.   0.21 -0.48
 1.01 -1.87  1.01  0.    0.    0.    0.   0.39 -1.78
 0.    1.01 -1.87  1.01  0.    0.    0.   0.48 -2.69
 0.    0.    1.01 -1.87  1.01  0.    0.   0.55 -3.49
 0.    0.    0.    1.01 -1.87  1.01  0.   0.48 -2.69
 0.    0.    0.    0.    1.01 -1.87  1.01 0.39 -1.78
 0.    0.    0.    0.    0.    1.01 -1.88 0.21 -0.48

```

The output has been formatted for better readability. It shows the full $`7\times 9`$ Jacobian evaluated at the argument:

```
      x(1) = 1.72
      x(2) = 3.45
      x(3) = 4.16
      x(4) = 4.87
      x(5) = 4.16
      x(6) = 3.45
      x(7) = 1.72

      prm(1) = 1.3
      prm(2) = 0.245828

```

#### 3.1.1 Compressed Jacobian – Seeding

The full Jacobian computed in the preceding section is sparse, and its sparsity pattern can be visualized as follows.

$$\begin{array}{ccccccccc}& & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \\ & & & & & & & & \end{array}.$$

Here, $``$ stands for a nonzero entry, and blanks represent structural zero entries in the Jacobian. In other words, no dependence exists between the corresponding dependent and independent variables.

Curtis-Powell-Reid (CPR)  seeding is based on the idea that certain columns of the Jacobian can be merged to share storage. For example, column 1 and column 7 could share one column, thus resulting in a compressed version of the Jacobian. This implies that the sparsity pattern must be known in advance in order to exploit matrix compression techniques. Recall that the derivative code generated by ADIFOR always loops over the derivative components of all active variables.An active variable w can be characterized as follows. At some point in the program (1) the value of w depends on the value of some independent variables and (2) there is some dependent variable whose value depends on w. Variables that are not active are referred to as passive. Many of them are equal to zero, leading to predictably trivial multiplications that one would like to avoid. Therefore, instead of computing the Jacobian as a Jacobian times identity matrix product, one could try to compute a compressed Jacobian using a seed matrix with fewer columns than the identity. Since the number of independent variables is often very large, the size of the seed matrix can become much smaller, leading to a decreased complexity of the forward vector mode Jacobian computation.

In CPR seeding, one considers the column incidence graph of the Jacobian to try to determine a minimal vertex coloring. Whenever two vertices share the same color, the corresponding columns can be stored in the same column of the compressed Jacobian. In our example, the column incidence graph has the following structure.

Unfortunately, since the vertex coloring problem is known to be NP-complete in general, the use of heuristics is essential. The coloring in the example graph has been found “by inspection.” Different colors are represented by different vertex shapes.

The number $`\nu `$ of different colors used determines the number of columns in the CPR seed matrix. Its rows are Cartesian basis vectors in $`IR^\nu `$. Whenever two vertices share the same color, the corresponding rows in the seed matrix contain the same Cartesian basis vector. This leads to the following seed matrix for our example.

$$\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}$$

The “compressed:” section of the makefile in Appendix B.2 builds an executable that produces output similar to the following.

```
-1.88  1.01  0.    0.21 -0.48
 1.01 -1.87  1.01  0.39 -1.78
 1.01  1.01 -1.87  0.48 -2.69
-1.87  1.01  1.01  0.55 -3.49
 1.01 -1.87  1.01  0.48 -2.69
 1.01  1.01 -1.87  0.39 -1.78
-1.88  0.    1.01  0.21 -0.48

```

It shows the compressed $`7\times 5`$ Jacobian evaluated at the current argument. The reconstruction of the original Jacobian is a simple substitution process as described in \[9, Chapter 7\].

#### 3.1.2 Sparse Forward Mode – SparsLinC

We compute the Jacobian of f with respect to x. The sparsity pattern of the Jacobian need not be known a priori. Derivative components are sparse vectors of (index, value) pairs. The computational overhead of sparse vector arithmetic results from the index calculations. Structural sparsity of the extended Jacobian \[9, Chapter 2\] is exploited by avoiding trivial multiplications by zero. The decision about when to apply runtime sparsity methods depends on the problem structure.

The use of SparsLinC with our example can be described by the following steps:

1. Use the file explosion.cmp without change.
2. Add the entry AD\_FLAVOR=sparse to explosion.adf.
3. Write driver program as shown in Appendix B.5.
4. Generate g\_explosion.f by calling Adifor2.1 AD\_SCRIPT=explosion.adf, and copy it from output\_files to the current directory.
5. Compile g\_explosion.f and the driver program.
6. Link with ReqADIntrinsics-Linux86.o, libADIntrinsics-Linux86.a, and
libSparsLinC-Linux86.a.
7. Run the executable.

Let us have a closer look at the driver program. The differentiated subroutine

g\_expl(dim,parmax,x,g\_x,prm,f,g\_f)

is generated in g\_explosion.f, where g\_x and f\_x are integer arrays of dimension dim containing pointers to the corresponding sparse derivative objects. Seeding g\_x is performed in the driver by

```
      do 10 i=1,dim
        g_x(i)=0
        CALL DSPSD(g_x(i),i,1.d0,1)
 10   continue

```

Notice that g\_x(i) must be initialized properly before calling DSPSD. The call initializes the sparse derivative object pointed at by g\_x(i) as (i,1.d0). Consequently, after executing this loop, g\_x contains the sparse identity in $`IR^{\mathrm{𝚍𝚒𝚖}}.`$

The sparse derivative components are extracted by

```
      do 30 i=1,dim
        g_f(i)=0
        CALL DSPXSQ(indvec,valvec,dim,g_f(i),outlen,info)
 30   continue

```

where indvec is the index vector and valvec the corresponding value vector making up the sparse representation of the derivative object pointed at by g\_f(i). The parameter outlen is the number of nonzero entries in the ith row of the Jacobian. Consequently, the outlen first entries of indvec and valvec define the nonzero entries in the ith row of the Jacobian as index value pairs (indvec(j), valvec(j)).

The “sparse:” section of the makefile in Appendix B.2 builds an executable that produces output similar to the following.

```
(1, -1.88) (2,  1.01)
(1,  1.01) (2, -1.87) (3,  1.01)
(2,  1.01) (3, -1.87) (4,  1.01)
(3,  1.01) (4, -1.87) (4,  1.01)
(4,  1.01) (5, -1.87) (6,  1.01)
(5,  1.01) (6, -1.87) (7,  1.01)
(6,  1.01) (7, -1.88)

```

The sparse Jacobian is given by sparse row vectors whose nonzero entries are represented by (index, value) pairs. The reconstruction of the full Jacobian is straightforward. Refer to  for further information on how to use SparsLinC.

### 3.2 ADIC 1.1

All files involved in using ADIC with the C version of our example problem are shown in Appendix C. The Jacobian of F with respect to x is computed in forward vector mode. ADIC uses the derived data type DERIV\_TYPE to associate derivative components with active variables. The data type InactiveDouble is used to indicate that some floating-point variable is not active. The call of ad\_AD\_Init(dim) causes the derivative code to use only the first dim elements of the derivative components. The latter are vectors of length GRAD\_MAX (see Appendix C). The vector x is declared to contain the dim independent variables by calling

```
      ad_AD_SetIndepArray(x,dim);
      ad_AD_SetIndepDone();

```

The function value components of DERIV\_TYPE variables are accessed by means of the macro DERIV\_val.

ADIC generates a differentiated version of the subroutine explosion and names it ad\_explosion. The argument list remains unchanged because all AD-related information is encapsulated in the new data type DERIV\_TYPE. If the function returns a double, however, the differentiated function becomes a procedure, and the returned value becomes the first argument. The derivative components of scalar variables of this type can be extracted by calling the routine ad\_AD\_ExtractGrad. The actual function value can be extracted with ad\_AD\_ExtractVal. In the example we call ad\_AD\_ExtractGrad(jac,F\[i\]) inside a loop over the dim elements of the vector of dependent variables F. The auxiliary variable jac is declared as a passive vector of size dim and contains the ith row of the Jacobian.

A makefile is provided to build the executable explosion.ad in order to compute the first seven columns of the Jacobian shown in Section 3.1.

### 3.3 ADOL-C

The AD-tool ADOL-C is based on operator overloading. Using this technique, one can log for each operation during the program execution the operator and the variables that are involved. Hence, one obtains a new internal representation of the function evaluation. Based on the generated execution log, ADOL-C computes the desired derivatives.

To apply ADOL-C for derivative calculations, one first has to modify the evaluation program to record the internal representation called tape. This modification starts with including header-file(s) that introduce the new data types and functions. Here, the easiest way is to simply include adolc.h. Second, one has to define the part of the program for which one wants to compute the derivatives. From now on, this part is called the active section. The statement trace\_on(tag,keep); determines the beginning of the active section, the statement trace\_off(file); the end of the active section. The parameter tag identifies the function to be differentiated. Hence, several function representations can be kept at the same time. Because of the shortness of this introduction, the optional parameters keep and file are not explained here but are explained in the documentation . Third, one has to change the types of the independents to adoubles and mark them as independents using the overloaded operator $`<<=`$. Similarly, one must mark the dependents using the overloaded operator $`>>=`$. Finally, all variables that lie on the way from the independents to the dependents must be declared as adoubles. This step includes the generation of a new function explosion\_ad, which contains the same source code as before but in the interface the double-variables have to be changed to adouble-variables. The required modification of the original source code is illustrated by Appendix D.

#### 3.3.1 Jacobian Calculation with the Easy-to-Use Drivers

After the tapes are generated during the execution of the active section, the required derivative objects can be calculated. For that purpose ADOL-C provides a variety of easy-to-use drivers: gradient(..), jacobian(..), hessian(..), vec\_jac(..) computing the product vector times Jacobian,
jac\_vec(..) computing the product Jacobian times vector, and so on. Moreover, ADOL-C supplies routines evaluating the Taylor coefficient vectors and their Jacobians with respect to the current state vector of solution curves defined by ordinary differential equations. Furthermore, there are drivers for derivative tensors and for the differentiation of implicit and inverse functions.

If one wants to compute the Jacobian of $`F`$ with respect to $`x`$ and the parameters $`s`$ and $`t`$, one has to declare a variable storing the derivative information

double \**J = myalloc(dim,dim+parmax);

where myalloc is provided by ADOL-C to allocate a two-dimensional array. Then, the statement

jacobian(tag,dim,dim+parmax,v,J);

causes the computation of the Jacobian of the function representation contained in the tape with the number tag at the point v. For a consistency check, the second and third parameter determine the number of dependents and independents, respectively. Hence, only the two statements given above have to be inserted after trace\_off() in order to compute the Jacobian of $`F`$.

#### 3.3.2 Calculation of Sparsity Pattern

ADOL-C provides for the computation of sparsity patterns the driver

jac\_pat(tag,dim,dim+parmax,v,rb,cb,Jsp,option);

If one sets rb and cb to NULL, jac\_pat computes the sparsity structure of the complete Jacobian at the point where the tape was generated and stores it in the unsigned int-array Jsp. The corresponding statements that have to be added to the source code are shown in Appendix D. If a certain block structure of the Jacobian is known, the unsigned int-vectors rb and cb can be applied to describe a compressed form of the independent and dependent variables.

#### 3.3.3 Forward and Reverse Mode Using Low-Level Routines

Arbitrarily high-order derivatives can be calculated by using the low-level functions of ADOL-C for the forward and reverse mode of AD. These routines are explained in detail in a short reference available from the ADOL-C Web page. In this article, only the computation of first-order derivatives is sketched. If one wants the full Jacobian, it is preferable to use vector modes of AD. For that purpose, ADOL-C provides the drivers fov\_forward(…) and fov\_reverse(…). Here, fov stands for first-order vector. The other acronyms of the low-level routines have a corresponding meaning. The statement

fov\_forward(tag,dim,dim+parmax,p,v,X,Fp,Y);

computes the derivative object Y $`=F^{}(v)`$X for X $`IR^{\mathrm{𝖽𝗂𝗆}+\mathrm{𝗉𝖺𝗋𝗆𝖺𝗑}\times p}`$. The statement

fov\_reverse(tag,dim,dim+parmax,q,U,Z);

computes the derivative object Z $`=`$ U$`F^{}(v)`$ for U $`IR^{q\times \mathrm{𝖽𝗂𝗆}}`$. To prepare this reverse sweep, one has to call an appropriate forward routines, the choice of which is described in detail by the short reference mentioned above. The code segments required for computing the full Jacobian with the low-level routines are contained in Appendix D.

### 3.4 TAPENADE

Our last example uses the alpha version of TAPENADE 1.0 to illustrate the use of reverse-mode AD for computing the Jacobian of f with respect to x and prm. All files needed to run the example are described in Appendix E. The makefile in Appendix E.1 can be used to generate the executable explosion.ad. Running the latter results in the computation of the full Jacobian as in Section 3.1 but this time using the reverse mode of AD.

The following command-line parameters are used to call tapenade:

* -head”: The name of the head routine
* -vars”: The names of the independent variables
* -cl”: Switch indicating that reverse, or cotangent linear, mode is used

It generates derivative code in explcl.f, which must be compiled together with the driver program shown in Appendix E.2 and linked with the routines for storing and restoring the values of the tape as described in \[9, Chapter 2\].

The current version of TAPENADE does not provide the reverse vector mode. Hence, the Jacobian must be computed as a sequence of Jacobian transposed times vector products. This computation is done in the driver program by successively initializing the derivative components g\_f of f as the Cartesian basis vectors in $`IR^{\mathrm{𝚍𝚒𝚖}}.`$ Thus, the Jacobian is accumulated row by row at a computational complexity proportional to the number of dependent variables. This feature becomes particularly interesting in the case of large, single gradients. Refer to  for further details on reverse-mode AD.

## Appendix A Original Code

```
      subroutine bratu(dim,parmax,x,prm,F)
      integer    dim    , parmax
C     independent variables
      double precision x(dim), prm(parmax)
C     dependent variables
      double precision F(dim)
C
      integer i
      double precision h

      h = 2.0/(dim+1)
      F(1) = -2*x(1)+h*h*prm(1)/12.0*(1+10*exp(x(1)/(1.0+prm(2)*x(1))))
      F(2) = x(1)+h*h*prm(1)/12.0*exp(x(1)/(1.0+prm(2)*x(1)))

      do 1 i=2,dim-1
        F(i-1) = F(i-1)+x(i)+h*h*prm(1)/12.0*exp(x(i)/(1.0+prm(2)*x(i)))
        F(i) = F(i)-2*x(i)+h*h*prm(1)/1.2*exp(x(i)/(1.0+prm(2)*x(i)))
        F(i+1) = x(i)+h*h*prm(1)/12.0*exp(x(i)/(1.0+prm(2)*x(i)))
    1 continue

      F(dim-1) = F(dim-1)+x(dim)+h*h*prm(1)/12.0*exp(x(dim)/(1.0
     *           +prm(2)*x(dim)))
      F(dim) = F(dim)-2*x(dim)
      F(dim) = F(dim)+h*h*prm(1)/12.0*(1+10*exp(x(dim)/(1.0
     *         +prm(2)*x(dim))))
      end

```

## Appendix B ADIFOR 2.0

The following files can be downloaded from

http://www-unix.mcs.anl.gov/̃naumann/ad\_tools.html

```
explosion.adf
explosion.cmp
explosion.driver.compressed.f
explosion.driver.f
explosion.driver.sparse.f
explosion.f
explosion.sparse.adf
makefile

```

### B.1 explosion.cmp

The composition file lists the names of all files containing subroutines subject to differentiation. In our simple example there is just explosion.f.

### B.2 makefile

The makefile can be used for computing the full, compressed, and sparse Jacobians using ADIFOR 2.0. It also ensures a proper cleanup of all files that are generated automatically during this process.

```
AD_LIB=/home/uwe/ADTOOLS/ADIFOR2/ADIFOR2.0D.lib

dense:Ψ
ΨAdifor2.1 AD_SCRIPT=explosion.adf
Ψcp output_files/g_explosion.f .
Ψg77 -g -c g_explosion.f explosion.driver.f
Ψg77 -g -o explosion.ad.dense -L$(AD_LIB)/lib *.o \
Ψ$(AD_LIB)/lib/ReqADIntrinsics-Linux86.o \
Ψ-lADIntrinsics-Linux86

compressed:Ψ
ΨAdifor2.1 AD_SCRIPT=explosion.adf
Ψcp output_files/g_explosion.f .
Ψg77 -g -c g_explosion.f explosion.driver.compressed.f
Ψg77 -g -o explosion.ad.compressed -L$(AD_LIB)/lib *.o \
Ψ$(AD_LIB)/lib/ReqADIntrinsics-Linux86.o \
Ψ-lADIntrinsics-Linux86

sparse:
ΨAdifor2.1 AD_SCRIPT=explosion.sparse.adf
Ψcp output_files/g_explosion.f .
Ψg77 -g -c g_explosion.f explosion.driver.sparse.f
Ψg77 -g -o explosion.ad.sparse -L$(AD_LIB)/lib *.o \
Ψ$(AD_LIB)/lib/ReqADIntrinsics-Linux86.o \
Ψ-lADIntrinsics-Linux86 \
Ψ-lSparsLinC-Linux86

clean:Ψ
Ψrm -fr output_files
Ψrm -fr AD_cache
Ψrm *.o
Ψrm g_*
Ψrm explosion.ad.*
Ψrm *\~

```

### B.3 Jacobian

#### B.3.1 explosion.adf

The Jacobian of the output variable f with respect to the two input variables x and prm of the top-level routine expl is computed. The parameter AD\_PMAX must be set to the total number of independent variables, that is, 9<sup>§</sup><sup>§</sup>§This is because the Jacobian is computed by forward vector mode with a seed matrix that is equal to the identity in $`IR^9.`$.

```
AD_PROG = explosion.cmp
AD_TOP  = expl
AD_PMAX = 9
AD_IVARS= x prm
AD_DVARS= f

```

#### B.3.2 explosion.driver.f

```
      program main
      implicit none
C
C     Example: Explosion Equation
C     Driver for computing Jacobian
C
      integer dim, parmax, n
      parameter (dim=7, parmax=2, n=9)
      integer i,j
C     independent variables
      double precision x(dim), prm(parmax)
C     derivative components of independent variables
      double precision g_x(n,dim)
      double precision g_prm(n,parmax)
C     dependent variables
      double precision f(dim)
C     derivative components of dependent variables
      double precision g_f(n,dim)

C     Initialization of input variables
      x(1) = 1.72
      x(2) = 3.45
      x(3) = 4.16
      x(4) = 4.87
      x(5) = 4.16
      x(6) = 3.45
      x(7) = 1.72
      prm(1) = 1.3
      prm(2) = 0.245828

C     Seeding (identity)
      do 20 i=1,n
        do 10 j=1,parmax
          if (i.eq.j+dim) then
            g_prm(i,j)=1.0
          else
            g_prm(i,j)=0.0
          endif
 10     continue
 20   continue

      do 40 i=1,n
        do 30 j=1,dim
          if (i.eq.j) then
            g_x(i,j)=1.0
          else
            g_x(i,j)=0.0
          endif
 30     continue
 40   continue

C     call differentiated subroutine
      call g_expl(n,dim,parmax,x,g_x,n,prm,g_prm,n,f,g_f,n);

C     print Jacobian
      do 60 i=1,dim
        do 50 j=1,n
          print*, "f’(", i, ",", j, ")=",g_f(j,i)
 50     continue
 60   continue

      end

```

### B.4 Compressed Jacobian

#### B.4.1 explosion.adf

This is the same as in the dense case.

#### B.4.2 explosion.driver.compressed.f

```
      program main
      implicit none
C
C     Example: Explosion Equation
C     Driver for computing compressed Jacobian
C
      integer dim, parmax, n
      parameter (dim=7, parmax=2, n=5)
      integer i,j
C     independent variables
      double precision x(dim), prm(parmax)
C     derivative components of independent variables
      double precision g_x(n,dim)
      double precision g_prm(n,parmax)
C     dependent variables
      double precision f(dim)
C     derivative components of dependent variables
      double precision g_f(n,dim)

C     Initialization of input variables
      x(1) = 1.72
      x(2) = 3.45
      x(3) = 4.16
      x(4) = 4.87
      x(5) = 4.16
      x(6) = 3.45
      x(7) = 1.72
      prm(1) = 1.3
      prm(2) = 0.245828

C     Seeding (CPR)
      do 20 i=1,n
        do 10 j=1,parmax
          g_prm(i,j)=0.0
 10     continue
 20   continue

      g_prm(n-1,1)=1.0
      g_prm(n,2)=1.0

      do 40 i=1,n
        do 30 j=1,dim
          g_x(i,j)=0.0
 30     continue
 40   continue

      g_x(1,1)=1.0
      g_x(1,4)=1.0
      g_x(1,7)=1.0
      g_x(2,2)=1.0
      g_x(2,5)=1.0
      g_x(3,3)=1.0
      g_x(3,6)=1.0

C     call differentiated subroutine
      call g_expl(n,dim,parmax,x,g_x,n,prm,g_prm,n,f,g_f,n);

C     print compressed Jacobian
      do 60 i=1,dim
        do 50 j=1,n
          print*, "f’(", i, ",", j, ")=",g_f(j,i)
 50     continue
 60   continue

      end

```

### B.5 Sparse Jacobian

#### B.5.1 explosion.sparse.adf

In order to make ADIFOR generate derivative code that can use the sparse forward mode provided by SparsLinC, the parameter AD\_FLAVOR must be set to sparse.

```
AD_PROG = explosion.cmp
AD_FLAVOR = sparse
AD_TOP  = expl
AD_PMAX = 9
AD_IVARS= x
AD_DVARS= f

```

#### B.5.2 explosion.driver.sparse.f

```
      program main
      implicit none
C
C     Example: Explosion Equation
C     Driver for computing sparse Jacobian using SparsLinC
C
      integer dim, parmax, n
      parameter (dim=7, parmax=2, n=9)
      integer i,j
C     independent variables
      double precision x(dim)
C     passive inputs
      double precision prm(parmax)
C     pointers to sparse derivative components of
C     independent variables
      integer g_x(dim)
C     dependent variables
      double precision f(dim)
C     pointers to sparse derivative components of
C     dependent variables
      integer g_f(dim)

C     (index, value) pairs
      integer indexes(dim)
      double precision values(dim)

C     values used in extraction routine
      integer outlen, info

C     Initialization of input variables
      x(1) = 1.72
      x(2) = 3.45
      x(3) = 4.16
      x(4) = 4.87
      x(5) = 4.16
      x(6) = 3.45
      x(7) = 1.72
      prm(1) = 1.3
      prm(2) = 0.245828

C     Initialization of sparse data structures
      call XSPINI

C     Seeding (sparse identity)
      do 10 i=1,dim
        g_x(i)=0
        g_f(i)=0
        call DSPSD(g_x(i),i,1.d0,1)
 10   continue

C     Call differentiated subroutine

      call g_expl(dim,parmax,x,g_x,prm,f,g_f)

C     Extract derivative components
      do 30 i=1,dim
        call DSPXSQ(indexes,values,dim,g_f(i),outlen,info)
        if (info.eq.0) then
          do 20 j=1,outlen
            print*, "indexes(", i, ",", j, ")=", indexes(j)
            print*, "values(", i, ",", j, ")=", values(j)
 20       continue
        endif
 30   continue

      end

```

## Appendix C ADIC 1.1

The following files can be downloaded from

http://www-unix.mcs.anl.gov/̃naumann/ad\_tools.html

```
explosion.c
explosion.driver.c
explosion.init
makefile

```

### C.1 makefile

```
AD_INC = -I$(ADIC)/include -I.
AD_LIB = -L$(ADIC)/lib/$(ADIC_ARCH)

all:
ΨadiC -d gradient -i explosion.init
Ψg++ -g -o explosion.ad $(AD_INC) $(AD_LIB) \
Ψexplosion.ad.c explosion.driver.c \
Ψ-lADIntrinsics-C -laif_grad -lm

clean:
Ψrm explosion.ad*
Ψrm ad_deriv.h
Ψrm *\~

```

The environment variable ADIC must be set to the directory in which ADIC has been installed.

### C.2 explosion.init

```
[SOURCE_FILES]
  explosion.c

[gradient]
  GRAD_MAX=7

```

The script file specifies the input files (only one in this case) and a variety of other parameters that can be looked up in . The definition of GRAD\_MAX, the length of the derivative components, is important. Its default value is 5, which would cause trouble in our case where dim=$`7>5.`$ Refer to  for other ways to set the value of GRAD\_MAX.

### C.3 explosion.driver.c

```
#include "ad_deriv.h"
#include <stdio.h>
#include <iostream>
#include <stdlib.h>

extern void ad_explosion(int, DERIV_TYPE*, DERIV_TYPE*, DERIV_TYPE*);

void main()
{
  int i,j;
  int dim=7;
  int parmax=2;
  // independent variables
  DERIV_TYPE *x = new DERIV_TYPE[dim];
  DERIV_TYPE *prm = new DERIV_TYPE[parmax];
  // independent variables
  DERIV_TYPE *F = new DERIV_TYPE[dim];
  // independent variables
  InactiveDouble *jac = new InactiveDouble[dim];

  ad_AD_Init(dim);

  ad_AD_SetIndepArray(x,dim);
  ad_AD_SetIndepDone();

  DERIV_val(x[0]) = 1.72;
  DERIV_val(x[1]) = 3.45;
  DERIV_val(x[2]) = 4.16;
  DERIV_val(x[3]) = 4.87;
  DERIV_val(x[4]) = 4.16;
  DERIV_val(x[5]) = 3.45;
  DERIV_val(x[6]) = 1.72;
  DERIV_val(prm[0]) = 1.3;
  DERIV_val(prm[1]) = 0.245828;

  ad_explosion(dim,x,prm,F);

  for (i=0;i<dim;i++) {
    ad_AD_ExtractGrad(jac,F[i]);
    for (j=0;j<dim;j++)
      cout << "f’[" << i+1 << "," << j+1 << "]=" << jac[j] << endl;
  }

  ad_AD_Final();
}

```

## Appendix D ADOL-C

### D.1 Modifications of Original Source Code

```
#include <adolc.h>
#include <SPARSE/sparse.h>
#include <stdio.h>
#include <iostream>
#include <stdlib.h>

extern void explosion_ad(int, adouble*, adouble*, adouble*);

void main()
{
  int i,j;
  int dim=7;
  int parmax=2;
  int tag = 1;

  // independent variables (passiv)
  double *v = new double[dim+parmax];
  // dependent variables (passiv)
  double *Fp = new double[dim];

  // independent variables (active)
  adouble *x = new adouble[dim];
  adouble *prm = new adouble[parmax];
  // dependent variables (active)
  adouble *F = new adouble[dim];

  v[0] = 1.72;
  v[1] = 3.45;
  v[2] = 4.16;
  v[3] = 4.87;
  v[4] = 4.16;
  v[5] = 3.45;
  v[6] = 1.72;
  v[7] = 1.3;
  v[8] = 0.245828;

  trace_on(tag);
    for(j=0;j<dim;j++)
      x[j] <<= v[j];
    for(j=0;j<parmax;j++)
      prm[j] <<= v[j];

    explosion_ad(dim,x,prm,F);

    for(j=0;j<dim;j++)
      F[j] >>= Fp[j];
  trace_off();

}

```

### D.2 Computation of Sparsity Pattern

```
  // Sparsity pattern declarations
  int option[3];
  unsigned int** Jsp = new unsigned int*[dim];

  for(j=0;j<dim;j++)
    Jsp[j] = new unsigned int[dim+parmax];
  option[0] = 0; // automatic detection for AD mode
  option[1] = 0; // save propagation of bit-pattern
  option[2] = 0; // no output

  ...

  // Sparsety pattern computation after trace_off

  jac_pat(tag,dim,dim+parmax,v,NULL,NULL,Jsp,option);

```

### D.3 Jacobian Calculation Using Low-Level Routines

```
  // Calculate Jacobian using forward mode driver
  double** X = myalloc(dim+parmax,dim+parmax);
  // Calculate Jacobian using reverse mode driver
  double** U = myalloc(dim,dim);

  ...

  // Use low level routines to compute Jacobian
  // forward:

  for (j=0;j<dim+parmax;j++)
  {
   for (i=0;i<dim+parmax;i++)
     X[j][i] = 0.0;
   X[j][j] = 1.0;
  }

  // first order vector mode forward
  // ^     ^     ^           ^^^^^^^ = fov_forward
  fov_forward(tag,dim,dim+parmax,dim+parmax,v,X,Fp,J);

  // reverse
  for (j=0;j<dim;j++)
  {
   for (i=0;i<dim;i++)
     U[j][i] = 0.0;
   U[j][j] = 1.0;
  }

  // prepare reverse sweep with appropriate forward sweep
  zos_forward(tag,dim,dim+parmax,1,v,Fp);

  // first order vector mode reverse
  // ^     ^     ^           ^^^^^^^ = fov_reverse
  fov_reverse(tag,dim,dim+parmax,dim,U,J);

```

## Appendix E TAPENADE

The following files can be downloaded from

http://www-unix.mcs.anl.gov/̃naumann/ad\_tools.html

```
explosion.f
explosion.driver.f
makefile

```

### E.1 makefile

```
ADSTACK = $(HOME)/ADTOOLS/TAPENADE/tapenade1.0/stack/adStack.o

all:
Ψtapenade -head expl -vars "x prm" -cl explosion.f
Ψg77 -o explosion.ad explcl.f explosion.driver.f \
Ψ$(ADSTACK)

clean:
Ψrm explosion.ad
Ψrm explcl.f
Ψrm -fr diffgen
Ψrm *\~

```

### E.2 explosion.driver.f

```
      program main
      implicit none
C
C     Example: Explosion Equation
C
C     Driver for computing Jacobian using TAPENADE’s
C     Reverse Mode
C

      integer dim, parmax, n
      parameter (dim=7, parmax=2, n=9)
      integer i,j
C     independent variables
      double precision x(dim), prm(parmax)
C     derivative components of independent variables
      double precision g_x(dim)
      double precision g_prm(parmax)
C     dependent variables
      double precision f(dim)
C     derivative components of dependent variables
      double precision g_f(dim)
C     the whole Jacobian matrix
      double precision jac(dim,n)

C     Initialization of input variables
      x(1) = 1.72
      x(2) = 3.45
      x(3) = 4.16
      x(4) = 4.87
      x(5) = 4.16
      x(6) = 3.45
      x(7) = 1.72

      prm(1) = 1.3
      prm(2) = 0.245828

C     Compute Jacobian as sequence of Jacobian transposed
C     times vector products
      do 40 i=1,dim
        do 10 j=1,dim
          g_f(j)=0.d0
          g_x(j)=0.d0
 10     continue
        g_f(i)=1.d0
        g_prm(1)=0.d0
        g_prm(2)=0.d0
        call explcl(dim,parmax,x,g_x,prm,g_prm,f,g_f);
        do 20 j=1,dim
          jac(i,j)=g_x(j)
 20     continue
        do 30 j=1,parmax
          jac(i,j+dim)=g_prm(j)
 30     continue
 40   continue

C     print Jacobian
      do 60 i=1,dim
        do 50 j=1,n
          print*, "f’(", i, ",", j, ")=",jac(i,j)
 50     continue
 60   continue

      end

```


text:
# 𝑋𝑌 checkerboard antiferromagnet in external field

## 1 Introduction

A huge degeneracy of classical ground states in geometrically frustrated magnets can be lifted by quantum or thermal fluctuations via a so called order by disorder effect . In the present work we study the thermal order by disorder effect for the $`XY`$ antiferromagnet on a checkerboard lattice. This lattice is a two-dimensional network of corner-sharing squares with crossings, which are topologically equivalent to tetrahedra, see Fig. 1. The present model can, therefore, be relevant to real pyrochlores with the easy-plane type anisotropy Er2Ti2O7 and Er2Sn2O7  and also have an experimental realization as an array of the Josephson junctions or a superconducting wires network in transverse magnetic field .

The considered model is described by the Hamiltonian

$$\widehat{}=J\underset{ij}{}𝐒_i𝐒_j\underset{i}{}𝐇𝐒_i,$$
(1)
where $`𝐒_i=(\mathrm{cos}\phi _i,\mathrm{sin}\phi _i)`$ are the classical planar spins and $`𝐇=(H,0)`$ is in-plane magnetic field. The sum is over nearest-neighbors pairs on a checkerboard lattice and $`J>0`$ is an antiferromagnetic coupling constant. The Hamiltonian can be expressed as a sum over elementary plaquettes (squares with crossings):

$$\widehat{}=\frac{1}{2}\underset{p}{\overset{N_p}{}}\left(J𝐒_p^2𝐇𝐒_p\right)JN,$$
(2)
where $`𝐒_p=_{ip}𝐒_i`$ is the total magnetisation of a plaquette, $`N`$ is the number of sites on the lattice and $`N_p=\frac{1}{2}N`$ is the number of plaquettes. Minimizing energy of a single plaquette one obtains the classical constraint

$$𝐒_p=𝐇/(2J).$$
(3)
The minimal total energy is reached if the above constraint is satisfied on every plaquette. In the field range $`0HH_{\mathrm{sat}}=8J`$ the ground state of the model remains underconstrained and infinitely degenerate with a finite entropy. This huge degeneracy is determined by two factors: (i) possible basis spin quartets, which obey the classical constraint, are parametrized by two continuous variables and (ii) by various periodic and aperiodic spin structures formed with the same basis quartet. The latter property can be seen by looking at the Fourier transform of the interactions on the checkerboard lattice , which has a flat momentum-independent branch with minimal energy.

## 2 Thermal order by disorder effect

Since the degeneracy of the classical ground state is a consequence of lattice topology rather than being a symmetry imposed property, various classical ground state configurations have different excitation spectra. At finite temperatures magnetic system spans a phase volume in the vicinity of the ground state manifold. Due to a varying density of excitations, the system can be effectively trapped in the neighbourhood of certain ground states. Such an ergodicity breaking leads to a lifting of zero-$`T`$ degeneracy of a frustrated magnet and to thermal order by disorder selection. The statistical weights of different ground state configurations $`\psi `$ are given by $`w[\psi ]\mathrm{exp}(F[\psi ]/T)`$, where $`F[\psi ]`$ is the partial free energy obtained by integrating out ‘fast’ excitation modes. The minimum of $`F[\psi ]`$ ensures that a macroscopic system is trapped in the vicinity $`\psi ^{\mathrm{min}}`$ and probability to find it in another classical ground state is vanishingly small.

To find which spin configurations are favoured by thermal fluctuations we start with a simple perturbative treatment over the mean-field result. In the ground state configuration a magnitude of a local field derived from Eq. (2) is the same $`H_{\mathrm{loc}}=2J`$ for all sites and all external fields $`0HH_{\mathrm{sat}}`$. The harmonic spin-wave Hamiltonian is expressed as

$$\widehat{}_2=H_{\mathrm{loc}}\underset{i}{}S_i^x+J\underset{ij}{}S_i^yS_j^y\mathrm{cos}\theta _{ij},S_i^x1\frac{1}{2}S_i^{y2},$$
(4)
where components of every spin are taken in its local coordinate frame and $`\theta _{ij}`$ is an angle between neighbouring spins. The first term in $`\widehat{}_2`$, which describes uncorrelated fluctuations of individual spins with $`S_i^{y2}=T/H_{\mathrm{loc}}`$, is taken as an unperturbed spin-fluctuation Hamiltonian, whereas the second term is a perturbation $`\widehat{V}`$. The correction to the free-energy is given by $`\mathrm{\Delta }F=\widehat{V}^2/2T`$:

$$\mathrm{\Delta }F=(T/8)\underset{ij}{}\mathrm{cos}^2\theta _{ij}=(T/8)\underset{ij}{}\left(𝐒_i𝐒_j\right)^2.$$
(5)
Thermal fluctuations produce, therefore, an effective biquadratic exchange, which lifts the zero-$`T`$ degeneracy in favour of maximally collinear states with the largest number of $`\mathrm{cos}\theta _{ij}=\pm 1`$. The harmonic Hamiltonian can be, of course, diagonalized with the help of the Fourier transform to obtain the spin wave modes $`\omega _𝐤^n`$ and their contribution to the free-energy:

$$\mathrm{\Delta }F_2=T\underset{n}{}\underset{𝐤}{}\mathrm{ln}(\omega _𝐤^n/\pi T).$$
(6)
The problem is reduced to minimization of the sum in (6). Still the real-space second order perturbation result appears to give a correct first insight especially for multisublattice configurations with $`n>2`$, when exact diagonalization of $`\widehat{}_2`$ becomes increasingly cumbersome.

An additional complication of highly frustrated magnets stems from a presence of several branches of zero modes: $`\omega _𝐤^m0`$, in which case Eq. (6) is not anymore correct. Instead the low-$`T`$ contribution to the free energy becomes:

$$\mathrm{\Delta }F=(N_2/2+N_4/4)T\mathrm{ln}(J/T)+T\underset{nm,𝐤}{}\mathrm{ln}\omega _𝐤^n/J+Tf_4,$$
(7)
where $`N_2`$ is the number of usual harmonic or quadratic modes, $`N_4`$ is the number of zero or soft modes, and $`f_4`$ is a contribution from interaction between soft modes and their interaction with quadratic modes. In order to estimate the last term one has to solve a nonlinear problem, which is, generally, a very complicated task. However, a partial selection between various spin configuration can be done on a basis of the leading $`T\mathrm{ln}(J/T)`$ term. If the total number $`N_2+N_4`$ is fixed ($`=N`$ for the $`XY`$ checkerboard antiferromagnet), then, the free energy is minimal for states with the maximum number of soft modes. The number of soft modes can be found either from direct diagonalization of $`\widehat{}_2`$ or from geometric consideration, which assigns a local soft mode to every void (empty square) with all spins around it being parallel or antiparallel to each other . Thus, the soft modes act similar to harmonic excitations and stabilize collinear states. The presence of soft modes is most easily seen in the low-temperature specific heat, which is found from Eq. (7) to be $`C=\frac{1}{2}N_4/4N`$. The deviation of the specific heat from a universal value of $`C=\frac{1}{2}`$ tells how many soft modes exist in a low-temperature state of the $`XY`$ checkerboard antiferromagnet.

## 3 Zero field behaviour

At zero magnetic field the ground state constraint (3) specifies configurations with $`𝐒_p=0`$ on every plaquette. Such configurations can be constructed either from noncollinear or collinear spin quartets. As was argued above, thermal fluctuations tend to select maximally collinear states with two up- and two down-spins on an arbitrarily chosen axis in the $`XY`$ plane. The gauge transformation $`S_i^{y\mathrm{down}}=S_i^{y\mathrm{down}}`$maps $`\widehat{}_2`$ for an arbitrary collinear state on the same reference harmonic Hamiltonian. Every collinear state has, therefore, the same harmonic spectrum and the same number of zero (soft) modes $`N_4=\frac{1}{2}N`$. Hence, the specific heat of a collinear state is $`C=\frac{3}{8}`$. All terms in the low-$`T`$ expression for the free energy (7) except the last one coincide for all collinear states. Our estimate indicates that the Néel state with the ordering wave-vector $`𝐐=(\pi ,\pi )`$ on the original square lattice has the lowest anharmonic contribution $`f_4`$ among all translationally symmetric states. This, however, does not necessarily mean an appearance of a quasi long-range order at $`𝐪=𝐐`$ when $`T0`$. So called weathervane defects  or wondering (rough) domain walls  have been considered as a source of disorder for the kagome antiferromagnet. For the $`XY`$ checkerboard antiferromagnet there are no zero-width domain walls and only weathervane defects can destroy a quasi long-range translational order. They cost zero classical energy and increase the free-energy by $`\mathrm{\Delta }F_d\epsilon _dT`$. A finite $`T`$-independent density of such defects can be estimated as $`n_d(1+e^{\epsilon _d})^1`$ by neglecting interaction between defects. These defects can destroy the true long-range order if their concentration exceeds the percolation threshold on the corresponding lattice. In this case only nematic correlations desribed by a traceless second-rank tensor $`O_{\mathrm{nem}}^{\alpha \beta }=S_i^\alpha S_i^\beta \frac{1}{2}\delta ^{\alpha \beta }`$ will be present at low temperatures.

In the absence of analytical theory to deal with this sort of behaviour we have tried to derive the necessary information from Monte Carlo simulations. We have calculated squares of the two relevant low-$`T`$ order parameters:

$$S(𝐐)=\frac{1}{N^2}\underset{i,j}{}S_i^\alpha S_j^\alpha e^{i𝐐(𝐫_i𝐫_j)},O_{\mathrm{nem}}^2=\frac{1}{N^2}\underset{i,j}{}S_i^\alpha S_i^\beta S_j^\alpha S_j^\beta \frac{1}{2}.$$
(8)
The results are presented in Fig. 1. There is a clear signature of the nematic order seen by the enhancement of $`O_{\mathrm{nem}}^2(N)`$ at low temperatures. The Néel order parameter reaches only very small values and does not show any appreciable enhancement. Considering a change in the scaling behaviour of the nematic order parameter $`O_{\mathrm{nem}}^2(N)1/N^\alpha `$ from $`\alpha =1`$, short-range correlations, at high temperatures to $`\alpha <1`$, power-law correlations, at low temperatures we estimate the Kosterlits-Touless transition temperature as $`T_{\mathrm{KT}}=0.014(2)J`$.

## 4 Finite field phases

In external magnetic fields $`0<H<H_{\mathrm{sat}}`$ the ground states with intermediate magnetisation (3) are, generally, noncollinear, except for $`H=\frac{1}{2}H_{\mathrm{sat}}`$, when a collinear ‘$`uuud`$’ state, Fig. 2b, belongs to the ground state manifold. Selection of the $`uuud`$ states by thermal fluctuations leads to a 1/2-magnetization plateau, which is similar to a 1/3-plateau of a classical kagome antiferromagnet . At all other fields only a partial collinearity is possible in the classical ground state. Geometric consideration suggests two prime candidate states shown in Fig. 2. The first canted state, Fig. 2a, exists in the whole range of fields and does not break remaining spin reflection symmetry about the field direction. The second partially collinear state, Fig. 2c, with three identical sublattices appears only for $`H>\frac{1}{2}H_{\mathrm{sat}}`$ and does break the mirror symmetry. It is easy to check that the fluctuation induced biquadratic exchange (5) favours a ‘more’ collinear state with broken reflection symmetry $`\mathrm{\Delta }F=(T/8)[3+\frac{1}{3}(8h^25)^2]`$, $`h=H/H_{\mathrm{sat}}`$ over a ‘less’ collinear canted state $`\mathrm{\Delta }F=(T/8)[2+4(2h^21)^2]`$ in the whole range of existence of the former state.

The translational degeneracy of the high-field partially collinear states is similar to the degeneracy of the $`uuud`$ states and corresponds to the dimer coverings of a square lattice , the total number of such states being $`1.157^N`$. The problem of lifting translational degeneracy for the $`uuud`$ states is similar to the zero-field case discussed above. Also, the $`uuud`$ have the same specific heat $`C=\frac{3}{8}`$. For the partially collinear states there are additional features in the thermal order by disorder effect. A general partially collinear state does not have soft modes. Soft modes, corresponding to all parallel spins around an empty square, exist, nevertheless, for certain translational patterns. The maximum number of soft modes $`N_4=N/4`$ appears for a state shown in Fig. 2d, which corresponds to a columnar arrangement of effective dimers. Such a translational pattern has the lowest $`T\mathrm{ln}(J/T)`$ contribution to the free energy. However, for all temperatures $`T0.001J`$ accessible with our Monte Carlo code we did not find a nonvanishing value of the structure factor corresponding to the this state. Thus, there is no a conventional Ising order parameter $`S_i^y0`$, $`i`$ in the high-field partially collinear state. Instead an Ising reflection symmetry is broken by a unique triatic order parameter:

$$T_{yyy}=S_i^yS_i^yS_i^y.$$
(9)
For Monte Carlo simulations of the $`XY`$ chekcer-board antiferromagnet we have used the standard Metropolis algorithm discarding $`10^5`$ Monte Carlo steps per spin (MCS) to reach thermal equilibrium and, then, average observables over $`10^6`$$`10^7`$ MCS. The lattice sizes were up to $`N=1024`$. In order to estimate statistical errors, all results have been averaged over 10–20 runs. The field dependence of the specific heat $`C`$ is shown in Fig. 2 for $`T=0.002J`$. $`C`$ starts at value, which is very close to analytically predicted $`\frac{3}{8}`$. In applied field the system looses gradually soft modes up to $`H=0.5H_{\mathrm{sat}}`$, where $`uuud`$ state with $`C=\frac{3}{8}`$ appears. The peaks in the specific heat at $`H_{c1}=4.15J`$ and $`H_{c2}=6.8J`$ indicate the phase transitions to a high-field state with broken reflection symmetry. The order parameters for this state are presented in the bottom right panel of Fig. 2. The square of the triatic order parameter $`T_{yyy}`$ is nonzero between $`H_{c1}`$ and $`H_{c2}`$. It shows a good statistical averaging and very little finite size dependence: results for smaller clusters fall on top of the presented data for $`N=1024`$. The inset shows $`S^{yy}(q)`$ with $`𝐪=(\pi ,\pi )`$ on a new lattice built from one sort of tetrahedra, which should be nonzero if the long-range order corresponding to Fig. 2d is present. The obtained data show that $`S^{yy}(q)0`$ in the thermodynamic limit. Thus, the high-field state is described by the triatic order parameter $`T_{yyy}`$, which makes the checkerboard lattice to be similar to a classical kagome antiferromagnet . Our results suggest that the triatic state survives up to $`T0.015J`$.

In conclusion, we have investigated the low temperature phases of the $`XY`$ checkerboard antiferromagnet in strong external fields. We have found that thermal fluctuations stabilize interesting phases with new type of broken symmetries: nematic order at $`H=0`$ and triatic order (9) in the field range $`\frac{1}{2}H_{\mathrm{sat}}<H<H_{\mathrm{sat}}`$. Further theoretical investigations should focus on a unique $`H`$$`T`$ diagram of this model.

## Acknowledgments

The authors thank Theory Group of the Institute Laue-Langevin for usage of the computer facilities.

## References


text:
# 1 Introduction

## 1 Introduction

The quantum-mechanical understanding of the physics of actinide compounds presents a challenge due to the intricate nature of the partially filled $`5f`$-shell. Compared to the rare-earths, for which the $`4f`$-states are most often completely localized, e.g. exhibiting atomic-like multiplet structure, the $`5f`$ states in the actinides are less inert and can play a significant role in bonding, depending on the specific actinide element and the chemical environment. This is most convincingly demonstrated in the elemental metals, for which a localization transition occurs when going from Pu to Am. In the early actinides, Th, Pa, U, Np, and the $`\alpha `$-phase of Pu, the relatively delocalized 5$`f`$-electrons actively contribute to bonding, and their atomic volumes decrease in a parabolic fashion, similarly to the behavior seen across the transition metal series. In Am, the $`f`$-electron localization is accompanied by an abrupt $`50\%`$ increase in the atomic volume, and for the heavier elements, Cm, Bk and Cf, the specific volume either remains constant or decreases only slightly. Pu lies at the borderline, and its very complex phase diagram implies that the $`f`$-electron properties are of particularly intricate nature. Depending on the chemical properties of the ligands, the actinide compounds may exhibit different degrees of $`f`$-electron localization for the same actinide element.

Over the past 30 years, the local spin density (LSD) and semi-local (generalized gradient - GGA) approximations to density functional theory have proven very useful and accurate in describing bonding properties of solids with weakly correlated electrons, demonstrating that the cohesive energy data for the homogeneous electron gas, that underlie these approximations, are representative of the conduction states in real materials. However, when 4$`f`$-electrons are involved, the atomic picture with localized partially filled $`f`$-shells is usually a better starting point for calculations. The most well known extensions of LSD, capable of describing electron localization, include the self-interaction corrected (SIC)-LSD, LDA+U, and orbital polarization methods.

It is possible to get a reasonable description of rare-earth materials with the LSD method by including a partially occupied $`f`$-shell into the core and projecting out the $`f`$-degrees of freedom from the valence bands. In such calculations, a combination of density functional theory with input from experimental data is used to describe bonding electrons and the atomic $`f^n`$ configuration, respectively. The SIC-LSD method can be viewed as effectively including an integer number of $`f`$-electrons in the core, however without restricting the unoccupied $`f`$-degrees of freedom. The localized $`4f`$ electrons in the rare earth metals and compounds have been well described by the SIC-LSD method. Intermediate valent Yb compounds have been described as a localized $`f^{13}`$ configuration plus a narrow $`f`$-band state pinned to the Fermi level. The free Yb atom is divalent with a completely filled $`f^{14}`$ shell. Thus, the destabilization of the localized $`f`$-manifold, which occurs in the solid state, is described in the SIC-LSD method by introducing two kinds of $`f`$-electrons. An integer number of $`f`$-electrons are localized while a non-integer number of hybridized band-$`f`$ electrons is determined by the self-consistent position of the Fermi level. A similar picture has emerged from calculations applying the LDA+U approach to Tm compounds. The number of localized $`f`$-electrons leads to a definition of valency of the actinide ions, given as the integer number of electrons available for band formation. Due to a substantial $`f`$-character of the valence bands, this valency is not the same as that determined by the total $`f`$-electron count, which includes both localized and itinerant $`f`$-electrons, and which is usually non-integral. Therefore, the most stable actinide valency in a given compound is determined by the balance between the localization energy and band formation energy (hybridization energy).

In section 2 of the present paper, the SIC-LSD method is briefly described. In section 3, the results for selected cases are presented and discussed, notably the $`\delta `$-phase of Pu, the actinide monopnictides and monochalcogenides, the PuO2 compound, and the UX3 intermetallics. The paper is concluded in section 4.

## 2 The SIC-LSD scheme

The electronic configuration of the actinide atoms is \[Rn\]$`5f^36d^17s^2`$ for U,
\[Rn\]$`5f^46d^17s^2`$ for Np, \[Rn\]$`5f^67s^2`$ for Pu, \[Rn\]$`5f^77s^2`$ for Am, and \[Rn\]$`5f^76d^17s^2`$ for Cm. In the solid state, the relative proportions of $`s`$, $`d`$ and $`f`$ electrons will change due to hybridization and charge transfer. The $`f`$ orbitals participate in bonding through their overlap with the $`f`$ and $`d`$ orbitals on neighbouring actinide ions, as well as, with the valence orbitals of the ligands. Either of these interactions lead to a broadening of the atomic $`f`$-level into an $`f`$ resonance, which one might hope to describe in two opposite limits, either as a hybridized band (as in the standard LSD picture), or as an atomic-like delta-function (by including a partially filled $`f^n`$ shell into the atomic core and decoupling all the other $`f`$-degrees of freedom, i.e., completely ignoring a possible $`f`$-electron contribution to bonding). These two extremal scenarios are depicted schematically in Figs. 1a and 1b, respectively, while the SIC-LSD scenario, which can be viewed as an interpolation between the two, is displayed in Fig. 1c.

In the SIC-LSD approach  the LSD total energy functional is corrected for the spurious self-interaction of each occupied state $`\alpha `$:

$$E^{SIC}=E^{LSD}\underset{\alpha }{\overset{occ.}{}}\delta _\alpha ^{SIC},$$
(1)
where the self-interaction correction, $`\delta _\alpha ^{SIC}`$, for a given state $`\alpha `$, is defined as the sum of the Hartree and exchange-correlation energies:

$$\delta _\alpha ^{SIC}=U[n_\alpha ]+E_{xc}^{LSD}[n_\alpha ].$$
(2)
This correction vanishes for an itinerant state, and therefore the SIC-LSD functional for such a state coincides with the conventional LSD functional. To benefit from the self-interaction correction, an electron state needs to spatially localize, which costs band formation energy due to loss of hybridization. Whether this is favorable depends on the relative values of the hybridization energy and the self-interaction correction energy. Hence, the latter is identified with the localization energy. The rationale behind the functional in Eq. (1) is that for a delocalized electron the interaction with a given atom is well described by the mean-field LSD potential. In contrast, the appropriate potential for a localized electron, due to a large Wigner delay time, will be corrected for the fact that other electrons on that atom rearrange in response to the presence of this localized electron. The self-interaction correction depends on the spatial distribution of the $`f`$ orbital, while the hybridization energy depends on the overlap of a given $`f`$ orbital with the $`f`$ and $`d`$ orbitals on the neighbouring actinide sites, and valence orbitals on the ligand sites. The $`f`$-electron which has become localized, by benefiting from the self-interaction correction, can no longer hybridize with the conduction electron bands to give rise to any band-related features or valency fluctuations. However, the $`f`$-states which have not been explicitly localized can hybridize with the conduction electrons and form fully or partially occupied bands (see Fig. 1c).

In the SIC-LSD formulation one deals with two types of $`f`$ electrons, the localized and hybridized $`f`$ electrons, as first implicated by Gschneider in relation to rare earths. By assuming different $`f^n`$ configurations of localized electrons, various valency configurations can be realized and studied in detail. Within SIC-LSD the valency is defined as the integer number of actinide valence electrons which are available for band formation, i. e.:

$$N_{val}=ZN_{core}N_{SIC}.$$
(3)
Here $`Z`$ is the atomic number, $`N_{core}`$ is the number of core (and semi-core) electrons (which for actinides is 86), and $`N_{SIC}`$ is the number of localized $`f`$-electrons on the actinide sites. Thus, e.g. a trivalent configuration of the actinide ions U3+, Np3+, Pu3+, Am3+, and Cm3+, is realized by localizing three, ($`f^3`$ configuration), four ($`f^4`$), five ($`f^5`$), six ($`f^6`$), and seven ($`f^7`$) $`f`$ electrons on the respective actinide atoms. For a given $`f^n`$ configuration, the minimum in the total energy as a function of lattice parameter determines the equilibrium lattice constant. By comparing the total energy minima for different $`f^n`$ configurations, the global groundstate configuration and lattice constant can be determined. In selecting the $`f^n`$ configuration the Hund’s rules are usually followed by alligning spins and maximizing the orbital moment in the direction opposite to the spins (for less than half-filled shells, or in parallel to the spins for more than half-filled shells). During the iterations towards self-consistency the localized states are allowed to relax, although generally they do not change much.

The SIC-LSD scheme has been implemented within the tight-binding linear-muffin-tin orbitals (TB-LMTO) method. The actinide semi-core $`6s`$ and $`6p`$ states have been described with a separate energy panel. Spin-orbit coupling has been fully included in the self-consistency cycles. For simplicity, for systems discussed here, we have assumed ferromagnetic arrangement of the magnetic moments.

## 3 Results and Discussion

### 3.1 $`\delta `$-Plutonium

Electronic structure calculations treating $`f`$-electrons as band states describe quite succesfully the equilibrium volumes of the early actinide metals. In an early study of Am, Skriver et al. found the $`f`$-electron localization, signalled by the onset of spin-polarization, giving rise to an almost full, and hence nonbonding, spin polarized $`f^7`$ band. Also, the high pressure phases of Am have been succesfully described by the standard LSD theory.  Recently, the SIC-LSD method was applied to the series of actinide metals from Np to Fm, correctly describing the itinerant nature of Np, the trivalency of Am, Cm, Bk and Cf, and the shift to divalency in Es and Fm. Pu turns out to be the most delicate case, being situated on the borderline between the itinerant and well localized actinides. The groundstate $`\alpha `$-phase is well reproduced by LDA calculations, but the high temperature $`\delta `$-phase is peculiar. The crystal structure is high-symmetry fcc, it has the largest specific volume of all Pu allotropes (25 % larger than that of $`\alpha `$-Pu), and the thermal expansion coefficient is negative. It has long been recognized that these facts are connected to a higher degree of localization of the $`f`$-electrons in the $`\delta `$-phase, but on the other hand the volume is still $`16\%`$ smaller than that of Am.

We have investigated the $`\delta `$-phase of Pu within the SIC-LSD approach. The total energy of a number of localization scenarios are summarized in Figure 2. In constructing the localized $`f^n`$ shells we have considered either $`LS`$ coupling (black curves) or $`jj`$ coupling (dashed curves) of the $`f`$-states. In the former case, all $`f`$-states are taken to be eigenstates of $`l_z`$ and $`s_z`$, with $`s_z`$ quantum number $`1/2`$ corresponding to spin-up, in accordance with Hund’s first rule, and $`l_z`$ quantum numbers occupied in the sequence $`3,2,1,0,1,2,3`$ to comply with Hund’s second rule (for less than half filled shells). In the latter case, one-electron $`f`$-states are taken as eigenstates of $`j^2`$ and $`j_z`$, with $`j=5/2`$ and $`j_z`$ quantum number occupied in the sequence $`5/2,5/2,3/2,3/2,1/2,1/2`$. It is important to stress that these are only the starting configurations of the localized states. Ultimately, the localized states are determined self-consistently by the SIC-LSD total energy minimization, but in practice the symmetry of the initial state tends to be preserved during iterations towards self-consistency. In other words there are energy barriers for a SIC state to drastically alter its symmetry.

One notices a significant energy gain for the configurations with large spin localized $`f^n`$ shells. The lowest energy is found for a localized $`LS`$ coupled $`f^5`$ shell, corresponding to trivalent Pu atoms. The equilibrium volume is 218 $`a_0^3`$, which is 30 % larger than the experimental volume of $`\delta `$-Pu. Clearly, this is not the appropriate representative of $`\delta `$-Pu. One more localized $`f`$-electron leads to an even larger volume and also a larger total energy, while fewer localized electrons do lead to smaller equilibrium volume but also larger total energy. The best agreement with the experimental volume within the $`LS`$ coupling scheme is obtained for the $`f^2`$ localized scenario. The $`jj`$ coupling scheme leads to a completely different picture. In this case the scenarios with localized $`f^2`$, $`f^3`$, and $`f^4`$ shells are almost degenerate in energy, with $`f^3`$ having the lowest energy, and an equilibrium volume $`12\%`$ smaller than the experimental $`\delta `$-Pu volume. The conclusion to be drawn here is that the LSD provides a poor account of the energetics of $`\delta `$-Pu: the exchange energy gained by the formation of large aligned spins is overestimated and leads to the wrong representation of the ground state. By taking $`jj`$ coupled localized shells, one artificially turns the spin-density contribution to the total energy off, and an improved description is obtained. This does not mean that the $`jj`$ coupled Pu $`f^3`$ shell is the correct ground state of Pu. Rather, the study demonstrates that more complicated ground states are called for. Within the restricted one-electron picture the $`jj`$ coupled localized shell is a better representation of the true ground state. Firstly, the true ground state must describe appropriately the spin fluctuations leading to the quenching of the Pu moment, secondly it is also likely that fluctuations in the number of localized $`f`$-electrons are needed. It is important to stress that $`\delta `$-Pu is a special case in the actinide series. When going to Am, the $`LS`$ ground state obtained with the SIC-LSD approach is quite adequate, leading to a localized $`f^6`$, $`M_S=3`$, $`M_L=3`$, i.e. $`J=0`$ ground state. The equilibrium volume is $`8\%`$ larger than the experimental volume, which is acceptable. The $`jj`$ coupled ground state is also $`f^6`$, $`J=0`$, in this case with a volume only $`3\%`$ larger than the experimental volume, but there is not such a drastic difference between the two representations for the well localized Am case as for Pu. Hence, the failure of the SIC-LSD in describing the highly correlated $`\delta `$-phase of Pu has been traced back to the large magnetic moment on Pu, persisting in the SIC-LSD description. By artificially setting the exchange interaction to zero, a much improved lattice constant has been obtained, as also found by Refs.  and . Since experiments find Pu to be non-magnetic, one must conclude that the mean-field approaches of LSD and SIC-LSD overestimate the tendency towards magnetic moment formation, by not taking into account quantum fluctuations in the $`f`$-shell. Recently, Savrasov et al. have presented a promising way of treating dynamical fluctuations and applied it successfully to $`\delta `$-Pu.

### 3.2 Actinide Monopnictides and Monochalcogenides

In the actinide monopnictides and monochalcogenides, which all crystallize in the NaCl structure at ambient conditions, the actinide-actinide separations are larger than in the elemental metals, and the tendency towards $`f`$-electron localization can already be observed from Np compounds onwards.

Here, we present the SIC-LSD electronic structures calculations of the monopnictides and monochalcogenides of U, Np, Pu, Am, and Cm. Figure 3 displays the calculated actinide ground state configurations through the series.

The calculations reveal clear trends towards more and more actively bonding $`f`$-electrons for a) lighter actinides, and b) lighter ligands. For the lighter actinides, the $`f`$-orbitals are more extended leading to larger overlaps with their nearest neighbours and smaller self-interaction corrections, both of these effects are favoring band formation. For the lighter ligands, in particular N and O, both the volume is decreased and ionicity is larger, the first of these effects leading to larger direct actinide-actinide overlap, and the latter effect favoring charge transfer.

The Cm compounds are the most localized systems, all exhibiting Cm in the trivalent $`f^7`$ configuration. The $`f^7`$ shell is so stable that variations of the ligand cannot disrupt its stability and scenarios with either one more or less localized $`f`$-electron have distinctly higher energies. Trivalency prevails in the Am compounds, but the stability of the $`f^7`$ shell causes the divalent Am state to be important in AmTe and AmPo. In the Pu compounds, the trivalent state also dominates, but for the lighter ligands $`f`$-electron delocalization sets in. In the Np compounds the tetravalent state dominates, while in the U compounds pentavalent states occur for the lighter ligands.

The densities of states of the actinide arsenides are shown in Fig. 4, with both trivalent and tetravalent actinide ions. In the trivalent case, the non-localized $`f`$-degrees of freedom give rise to narrow unoccupied bands above the Fermi level. In the tetravalent case the additional delocalized $`f`$-electron appears as an extra $`f`$-band. In Cm, this band appears far below the Fermi level, while in Am, Pu and Np this band lies just below the Fermi level. The band formation energy due to this extra band is sufficiently large in NpAs to outweigh the localization energy, and the tetravalent configuration becomes the ground state.

### 3.3 Plutonium Dioxide

PuO2 is the most favored compound for storage of Pu from nuclear waste. In the stoichiometric compound Pu is tetravalent with a localized $`f^4`$ shell, filled O $`p`$ bands and a large insulating gap. Recently, the chemical inertness of PuO2 has been questioned, in particular whether reactions with water could lead to further oxidation and the formation of PuO2+x. We investigated the PuO2+x system with the SIC-LSD approach by constructing a supercell with four PuO2 units together with an additional O atom in the interstitial region, thus forming a model of a PuO2.25 compound. The interesting conclusion to be drawn from the total energy as a function of $`f`$-localization (fig. 5) is that the nearest neighbour Pu atoms of the interstitial O transform to the pentavalent configuration by delocalizing one $`f`$-electron, which is donated to the extra O to form hybridized states, occuring in the gap-region of the pure PuO2 compound. Similarly, an O vacancy in PuO2 will lead to the formation of trivalent $`f^5`$ Pu ions in the vicinity of the vacancy. In effect, the localized $`f^n`$ shell of Pu acts as a reservoir for absorbing or releasing electrons to be accomodated by the chemical bonds of the O atoms. The lattice constant of the PuO2+x system almost does not vary with $`x`$ due to two opposing effects (Fig. 5). The added O per se leads to lattice expansion, but the additional bonding due to the formation of pentavalent Pu causes the lattice to contract.

### 3.4 UX3 compounds, X=Rh, Pd, Pt, Au

The sequence of Uranium intermetallics, UX3, X=Rh, Pd, Pt and Au are interesting due to their variation in metallic properties. This is well reproduced in the SIC-LSD approach. The U configuration changes from $`f^0`$ in URh3 to $`f^1`$ in UPt3 and to $`f^2`$ in UPd3 and UAu3. This is due to hybridization of the ligand $`d`$-band, which in URh3 is not completely occupied, with the U $`f`$-electrons to form hybridized bands rather than non-bonding localized states. In UPd3 and UAu3 the $`d`$-band is full, and the U $`f`$-electrons can not contribute further to the bonding, while UPt3 is the borderline case, where the $`d`$-band is full but sufficiently close to the Fermi level that the $`f`$-electrons can hybridize in, i.e., the U $`f`$-manifold is situated in between the fully delocalized and fully localized scenarios, in good accord with the observation of heavy fermion properties of this compound.

## 4 Conclusions

In summary, we have reviewed the electronic properties of a number of $`f`$-electron systems, as obtained within SIC-LSD approach. We have demonstrated that this approach is well suited to describe trends regarding lattice parameters and valencies of these systems. It works especially well for systems with well localized $`f`$-shells. It also indicates that the ground state of $`\delta `$-Plutonium is more complex than the SIC-LSD can describe and thus underlines the need of developing a dynamic generalization of the SIC-LSD approach.

## 5 Acknowledgments

This work has been partially funded by the Training and Mobility Network on ‘Electronic Structure Calculation of Materials Properties and Processes for Industry and Basic Sciences’ (contract:FMRX-CT98-0178) and by the Research Training Network on ‘Ab-initio Computation of Electronic Properties of f-electron Materials’.


text:
# 1 A numerical estimate of the energy gap Δ for a chain with 𝑁= 9. Phase boundary is obtained from the location of Δ(Γ)=0, and is shown in the inset.

Quantum Spin Glass Phase Boundary in $`\pm J`$

Transverse Field Ising Systems111 Presented at CMDays-03, Jadavpur Univ., Kolkata, Aug. 2003. Proc. in Ind. J. Phys. (to be published).

Arnab Das, Amit Dutta and Bikas K. Chakrabarti

Saha Institute of Nuclear Physics, 1/AF Bidhannagar, Kolkata - 700064, India

Physics Department, Indian Institute of Technology, Kanpur - 208016, India

Abstract: Here we study zero temperature quantum phase transition driven by the transverse field for random $`\pm J`$ Ising model on chain and square lattice. We present some analytical results for one dimension and some numerical results for very small square lattice under periodic boundary condition. The numerical results are obtained employing exact diagonalization technique following Lanczos method.

Keywords: Quantum phase transition, Spin glasses, Transverse Ising model, Diagonalization techniques.

PACS Nos.: 05.70.Fh, 42.50.L, 75.10.N

1. Introduction

The interest in the study of transverse Ising spin glass models was revived in early 1990 by the discovery of zero-temperature transition in dipolar Ising transverse field magnet LiHoxY1-xF4 . Proton glasses such as mixture of ferroelectric and anti-ferroelectric materials like Rb(1-x)(NH4)x(H2P)4  also provided earlier useful realizations of such quantum spin glasses.

These developments initiated extensive theoretical studies in quantum spin glass models. Ising model in transverse field has already been studied extensively in this context through analytical approaches using approximate renormalization techniques and real space renormalization group method, as well as using numerical methods like quantum Monte Carlo and exact diagonalization techniques . Fairly extensive studies on the quantum spin glass phases have been made so far for transverse field Edward-Anderson model and Sherrington-Kirkpatrick model with random exchange distributions  and to some extent on quantum Heisenberg spin glass models . In view of the rigorous developments in the study of two dimensional nearest neighbour Ising model with random $`\pm J`$ exchange interactions and the precise knowledge of location of the Nishimori line  in such classical spin glass model (driven by temperature), we consider here the quantum phase transition (at zero temperature) in the same $`\pm J`$ Ising model.

We have shown analytically that introduction of random $`J`$ impurities cannot affect the zero temperature phase transition in one dimensional system as they can be transformed away. We have also compared and verified the result numerically for the small system size considered. For two dimensional systems, we present some priliminary results obtained for a square lattice using exact diagonalization results for very small system sizes following Lanczos technique . Only the behaviours of configurationally averaged energy gap $`\mathrm{\Delta }=(E_1E_0)`$ between the first excited state and the ground state and the second order response function $`\chi =\left(^2E_0/\mathrm{\Gamma }^2\right)`$, equivalent to specific heat, have been studied here. The variations of $`\mathrm{\Delta }`$ and $`\chi `$ with respect to transverse field $`\mathrm{\Gamma }`$ have been obtained, and the phase boudary has been estimated from these results.

We work with a transverse Ising system, using only nearest neighbour interactions, whose Hamiltonian is given by

$$=\underset{i,j}{}J_{ij}S_i^zS_j^z\mathrm{\Gamma }\underset{i=1}{\overset{N}{}}S_i^x,$$
(1)
where the transverse field $`\mathrm{\Gamma }`$ is uniform through out the system and the nearest neighbour exchange constant $`J_{ij}`$’s are chosen randomly from the binary distribution

$$P(J_{ij})=p\delta (J_{ij}+J)+(1p)\delta (J_{ij}J).$$
(2)
Here $`J`$ is taken positive and $`p`$ is thus the concentration of anti-ferromagnetic $`J`$ bonds in the system.

2. Results in One Dimensional System

Here first we show analytically that in a one dimensional transverse Ising Hamiltonian with uniform $`J`$ and $`\mathrm{\Gamma }`$, if we replace some $`J`$ bonds by $`J`$ bonds randomly, then the resulting Hamiltonian can be gauge transformed back to one with uniform $`J`$, and hence the critical field remains unchanged with randomness concentration. Simillar result for one dimensional system with distributed $`J`$ had been obtained earlier .

Let us take the one dimensional random bond tranverse Ising Hamiltonian

$$=\underset{i}{}J_iS_i^zS_{i+1}^z\mathrm{\Gamma }S_i^x,$$
(3)
where the transverse field $`\mathrm{\Gamma }`$ is uniform throughout the system, and $`J_i`$’s are randomly chosen from the same distribution as given in $`(2)`$. Since the $`J_i`$’s have same magnitude $`J`$ all through, and their randomness is only in their sign, we may write $`J_i=J\mathrm{sgn}(J_i)`$ , and thus Hamiltonian (3) takes the form

$$=J\mathrm{sgn}(J_i)S_i^zS_{i+1}^z\mathrm{\Gamma }S_i^x.$$
(4)
Now let us define a new set of spin variables as below

$$\stackrel{~}{S}_i^z=S_i^z\underset{k=1}{\overset{i1}{}}\mathrm{sgn}(J_k)$$

$$\stackrel{~}{S}_i^x=S_i^x$$

$$\stackrel{~}{S}_i^y=S_i^y\underset{k=1}{\overset{i1}{}}\mathrm{sgn}(J_k).$$

It is easy to see that $`\stackrel{~}{S}`$ ’s satisfy the same commutation and anti-commutation relations as those of $`S`$’s and hence will exhibit exactly the same dynamical behaviour. Now,

$$\stackrel{~}{S}_i^z\stackrel{~}{S}_{i+1}^z=S_i^zS_{i+1}^z\left[\underset{k=1}{\overset{i1}{}}[\mathrm{sgn}(J_k)]^2\right]\mathrm{sgn}(J_i),$$

or,

$$\stackrel{~}{S}_i^z\stackrel{~}{S}_{i+1}^z=S_i^zS_{i+1}^z\mathrm{sgn}(J_i),$$

since $`[\mathrm{sgn}(J_k)]^2=1`$. Thus in terms of new spin variables, Hamiltonian (4) becomes

$$=J\underset{i}{}\stackrel{~}{S}_i^z\stackrel{~}{S}_{i+1}^z\mathrm{\Gamma }\underset{i}{}\stackrel{~}{S}_i^x.$$
(5)
The above Hamiltonian describes the same random system in terms of new variables, and yet, as one can see, it has in itself no randomness at all. One can use Jordon-Wigner transformation in terms of $`\stackrel{~}{S}`$’s and see that here also quantum phase transition occurs only at $`\mathrm{\Gamma }\mathrm{\Gamma }_c(=J)`$ as it occurs in a non-random Hamiltonian in $`S`$’s. In Fig. 1, we present some data computed for a chain of size $`N=9`$, which shows that the gap $`\mathrm{\Delta }`$ vanishes at $`\mathrm{\Gamma }_c1`$ (the field being scaled by $`J`$). These data for $`\mathrm{\Delta }=E_1E_0`$ is obtained from the computed average value of $`E_0`$ and $`E_1`$, each one averaged over about 10 configurations for $`p0`$. For infinite system, $`\mathrm{\Delta }`$ is a linear function of $`\mathrm{\Gamma }`$ for $`\mathrm{\Gamma }\mathrm{\Gamma }_c`$. In our case, linearity is observed at high values of $`\mathrm{\Gamma }`$, and $`\mathrm{\Gamma }_c`$ is determined by backward linear extrapolation from the linear region.

From the numerical data in Fig. 1, we see that there is a slight variation of $`\mathrm{\Gamma }_c`$ with $`p`$ ($`\mathrm{\Gamma }_c`$ varies between 0.9 and 1.0). This variation can be attributed to the very small size of the system. However, it may be noted that with even number of $`J`$ bonds in the chain, with periodic boundary condition, there is no problem of incommensuration and $`E_0(\mathrm{\Gamma })`$ or $`E_1(\mathrm{\Gamma })`$ become strictly identical for such values of $`p`$. Similarly, in every case of odd number of $`J`$ bonds in the chain, incommensuration problem always occurs for one spin only, rendering identical values (but different from the even $`J`$ case) for $`E_0(\mathrm{\Gamma })`$ and $`E_1(\mathrm{\Gamma })`$ in all such cases.

3. Results for Square Lattice

We consider now the same system (represented by Hamiltonian (1)) on a square lattice of size 3$`\times `$3 with periodic boundary condition. We again calculate the ground state and the first excited state energy $`E_0`$ and $`E_1`$ respectievely as functions of the transverse field $`\mathrm{\Gamma }`$, for different values of $`p`$. Each value of $`E_0`$ and $`E_1`$ is averaged over at least 10 configurations for each $`p0`$. Apart from $`\mathrm{\Delta }`$, we also calculate $`\chi =^2E_0/\mathrm{\Gamma }^2`$ and their variations with $`\mathrm{\Gamma }`$ as shown in Figs. 2 and 3 respectively.

Our results here are severely constrained by the system size. The value of pure ferromagnetic critical field $`\mathrm{\Gamma }_c(p=0)`$ is found here to be around 2.2, while the series results  or cluster algorithms  give the value to be around 3.0. This discrepancy is attributed to the smallness of our system size $`(N=3^2)`$. However the qualitative behaviour of the order-disorder phase boundary (between ferro/spin glass and para) seems to be reasonable: $`\mathrm{\Gamma }_c(p)`$ decreases with $`p`$ initially, and then increases again as $`p`$ approaches unity (pure anti-ferromagnet). The use of periodic boundary condition here (to avoid some numerical errors) also restricts the domain features and thereby affects our results. The absence of the knowledge of the ground-state wave function (and the correlation functions) in this method also forbids us to analyse the structure of the ordered phases.

Acknowledgement: BKC is grateful to Hidetoshi Nishimori for useful discussions.

References

 W. Wu, B. Ellmann, T.F. Rosenbaum, G. Appeli and D.H. Reich, Phys. Rev. Lett. 67 2076 (1991); W. Wu, D. Bitko, T.F. Rosenbaum and G. Appeli, Phys. Rev. Lett. 71 1919 (1993)

 R. Pric, B. Tadic and R. Blinc, Z. Phys. B 61 69 (1985)

 B. K. Chakrabarti, A. Dutta and P. Sen, Quantum Ising Phases and Transitions in Transverse Ising Models, Lecture Notes in Physics, m41, Springer, Heidelberg (1996)

 R. N. Bhatt, in Spin Glasses and Random Fields, Ed. A. P. Young, p. 225, World Sc., Singapore (1998)

 H. Nishimori, Prog. Theor. Phys. 66 1169 (1980); Statistical Physics of Spin Glasses & Information Processing : An Introduction, Oxford University Press, Oxford (2001); H. Nishimori and K. Nemoto, J. Phys. Soc. Jpn. 71 1198 (2002)

 J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, Text in Appl. Maths. 12, Springer-Verlag, New York, (1993 )

 B. McCoy, in Phase Transitions and Critcal Phenomena, vol. II, Eds. C. Domb and M. S. Green, Academic Press, London (1983); D. S. Fisher, Phys. Rev. B 50 3799 (1994)

 H. Rieger and N. Kawashima, Euro. Phys. J. B 9 233 (1999)