# Shafranov shift correction to the Furth-Yoshikawa scaling of tokamak adiabatic compression

T Nicolas, Hinrich Lütjens, O Sauter, X Garbet

► To cite this version:

T Nicolas, Hinrich Lütjens, O Sauter, X Garbet. Shafranov shift correction to the Furth-Yoshikawa scaling of tokamak adiabatic compression. Plasma Physics and Controlled Fusion, 2023, 65 (5), pp.054004. 10.1088/1361-6587/acc5ae hal-04248723

HAL Id: hal-04248723

https://hal.science/hal-04248723v1

Submitted on 20 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Shafranov shift correction to the Furth-Yoshikawa scaling of tokamak adiabatic compression

T. Nicolas¹, H. Lütjens¹, O. Sauter², X. Garbet³

E-mail: timothee.nicolas@polytechnique.edu

¹CPHT, Ecole Polytechnique, CNRS, 91128 Palaiseau, France

²Swiss Plasma Center – EPFL, Lausanne, Switzerland

³CEA, IRFM, F-13108 Saint-Paul-Lèz-Durance, France

**Abstract.** In 1970, Furth and Yoshikawa[1] introduced the scalings of adiabatic plasma compression. Basically, if the shape of the external plasma boundary and the aspect ratio are preserved during the compression, then the density, kinetic pressure, beta and current scale respectively as $n \sim C^3$, $p \sim C^5$, $B \sim C^2$, $\beta \sim C$, $I_t \sim C$, where C is the linear compression ratio, that is, the ratio between initial and final major radii. In this work, we show analytically, by expanding the Grad-Shafranov equation in terms of C, that the deviation to the Furth-Yoshikawa scaling is related to the Shafranov shift that arises when beta increases at large compression ratios. There is an obvious effect of the Shafranov shift because the axis is moved to a region with larger volume element, and an indirect effect, associated to the relation between flux and radius. The latter effect adds to the first, and is of the same order of magnitude. The result is that the pressure increases less than the $C^5$ scaling, which can have a significant impact on the fusion power achieved at maximum compression. The analytical results are backed up by equilibrium simulations carried out with the CHEASE code. Equilibria are obtained for different values of C, with conservation of the total fluxes, $q$ profile, and entropy of the plasma. The agreement of the theory and simulations is very good when the boundary of the plasma is circular and the aspect ratio small. When the aspect ratio is close to 1, and/or the boundary not circular, the analytical result gives the gist of the reduction of compression. Finally, a pressure anisotropy $(p_{\perp} - p_{\parallel})/p$ approximately equal to the increase in normalized Shafranov shift is predicted.
# 1. Introduction

Magnetized target fusion intends to achieve a fusion relevant triple product by compressing a magnetized plasma using an imploding metal liner. In this paper, we focus on the proposal of the General Fusion startup, where the liner is a liquid lead-lithium eutectic, and the magnetized plasma a tokamak plasma [2]. The approach is largely based on the LINUS concept [3, 4], but in a toroidal configuration. Although the confinement device is very distinct from a tokamak, the plasma can be called a tokamak plasma, because it is (i) axisymmetric and (ii) comprises both toroidal and poloidal fields, with a safety factor similar to that of a tokamak. Although the compression of a tokamak plasma by a conducting liner has never been attempted outside of General Fusion, the adiabatic compression of tokamak plasmas has been investigated both theoretically and experimentally in the seventies. The purpose was to provide an alternative way of bringing the tokamak plasma to fusion conditions, by rapidly changing the external coil parameters ensuring the equilibrium, in order to produce an adiabatic compression of the plasma[1, 5, 6, 7, 8]. Consequently, the basic scalings of pressure, magnetic field, current and plasma $β$ during compression were provided in 1970 by Furth and Yoshikawa [1]. The scaling is derived as a consequence of particle, entropy and flux conservation. Several cases were considered regarding a possible difference in compression of minor radius $a$ and major radius $R_0$. Here, we will restrict ourselves to the case of constant aspect ratio $ε = a/R_0$, indeed we find a deviation to the scaling even in this simple case. In that case, by defining $C$ as the ratio between final and initial major radii, the scalings are as follows, for the pressure, density, magnetic field, toroidal current, total and poloidal beta:

$$ 
\begin{aligned}
p &\sim C^5 \\
n &\sim C^3 \\
B &\sim C^2 \\
I_t &\sim C \\
β &\sim C \\
β_P &\sim C
\end{aligned} 
\qquad (1) $$

The purpose of this paper is to present the deviation with respect to this scaling for finite $β$ values. Strictly speaking, the equilibrium calculations we perform are a small subset of the more general theory of circular cross section tokamak equilibrium originating in the work of Zakharov and Shafranov (see in particular volumes 2 and 11 of the Reviews of Plasma Physics, refs. [9, 10]). In particular, the case of tokamak compression with constant total flux and entropy (adiabatic compression) was examined extensively by Green et al. [5], using analytical developments rooted in the aforementioned pioneering work of Shafranov and Zakharov. The deviation to the scaling of Furth and Yoshikawa is, in a sense, contained in reference [5]. However, there is a difficulty with the convention used by the authors, which prevents one from deriving the final result, as explained in Appendix C. As a matter of fact, the subsequent works, refs. [11, 12, 13, 14], which simulated the adiabatic compression, observed and commented a deviation to the scaling, but failed to explain it, and do not even cite reference [5]. For example, in ref. [14], the deviation is attributed to the change of shape of the plasma boundary only. Moreover, we have examined the wealth of literature citing
ref. [5], and found none that discuss the Furth-Yoshikawa scalings directly. The deviation effect can be large and have a strong impact on fusion power when $\beta$ reaches several tens of percent, as is envisioned in the General Fusion proposal [2]. For these reasons, we felt it was worth revisiting the problem of adiabatic compression of tokamak plasmas. Because of the issue discussed in Appendix C, it turned out easier, both from a technical point of view and in order to ease the readability of the paper, to rederive the effect from scratch. In order to make it analytically tractable, a special case is considered (perfectly conducting wall with constant aspect ratio and shape, perfectly adiabatic).

In technical terms, the question that is addressed is as follows. Given an initial safety factor (or current) and pressure distribution for an ideal tokamak plasma limited by a wall assumed to be perfectly conducting, how do the profiles evolve when the wall is moved so as to provoke a reduction of the plasma volume at constant shape and aspect ratio? Of course, in the actual experiments attempted by General Fusion, the assumptions of perfect conductivity for the liquid metal wall and of adiabaticity for the plasma are unrealistic, and an actual prediction for the behaviour of the plasma in their experiments would require to take into account the wall and plasma resistivity, as well as the evolution of the shape of the boundary. The point, here, is to show that even in the most ideal case, assuming no transport whatsoever, the situation is worse than the optimistic conclusions based solely on the Furth and Yoshikawa scaling (1), which neglects plasma profiles. That is, the axis pressure (and, as a result, the fusion power, which scales as $p^2$) increases significantly less than predicted by the scaling, as seen, for example, in figure 6. The real plasma, which has losses of magnetic flux, particles and energy, will perform below this ideal performance.

The result is easily expressed in terms of a single parameter, the Shafranov shift. When the volume of the plasma is reduced, the adiabaticity law leads to the increase of the pressure, following the scaling (1) at lowest order. At the same time, the plasma $\beta$ increases (linearly with the linear compression ratio). The reason for that increase is that the magnetic pressure, which is linked to the flux conservation, increases less rapidly than the pressure. The $\beta$ increase leads to an increase of the Shafranov shift, which takes the magnetic axis to larger major radii, with larger specific volume. This larger volume is responsible for one part of the reduced pressure increase, the other part being due to an effect of the modification of the relation between poloidal flux and magnetic surface minor radius.

The analytical result obtained with circular flux surfaces is tested with equilibrium calculations in the case of both circular large aspect ratio and shaped tight aspect ratio equilibria. It is shown how a simple adaptation of the analytical result to the shaped case matches the CHEASE results in the range of parameters tested so far.

The paper is organized as follows. In section 2, the adiabaticity laws are written in local form, and the basic effect of the Shafranov shift is described. In section 3, the full analytical result in the case of large aspect ratio and circular cross section is derived. In section 4, comparisons with CHEASE simulations [15] are performed and a discussion follows in section 5 to conclude.
## 2. Physics of the compression

### 2.1. Adiabatic constraints

The plasma is supposed to undergo compression by the physical motion of a boundary made of a highly conducting material (eutectic Li-Pb liquid metal in the proposal of General Fusion [2]). If the compression is fast compared to the time scales of magnetic flux, particle and energy losses, then the compression is adiabatic. The time scales mentioned here pertain not only to the global loss time scales, but also to the transport time scales within the plasma. Assuming the plasma remains axisymmetric, the topology of the magnetic surfaces is preserved and the position in the plasma can at any time be described with the generalized coordinates $(\psi, \theta, \varphi)$, where $\theta$ and $\varphi$ are the poloidal and toroidal angles, and $\psi$ the poloidal flux divided by $2\pi$. The density and temperature are assumed to be flux surface quantities. Any surface $\psi = C^{\text{te}}$ contains at any time the same quantity of particles, entropy and magnetic flux (poloidal and toroidal). This leads immediately to the fact that the safety factor $q(\psi) = d\psi_T/d\psi$ is preserved, where $\psi_T$ is the toroidal flux divided by $2\pi$. The total number of particles contained within a surface bounded by $\psi$ is $N(\psi) = 2\pi \int d\psi n(\psi) \oint \mathcal{J} d\theta$, where $\mathcal{J} = [(\nabla\varphi \times \nabla\psi) \cdot \nabla\theta]^{-1}$ is the jacobian of the change of variables from cartesian to $(\psi, \theta, \varphi)$. The volume within a surface bounded by $\psi$ is $V(\psi) = 2\pi \int d\psi n(\psi) \oint \mathcal{J} d\theta$. The latter is often written $2\pi \int d\psi dl_p/B_p$, where $dl_p$ is the poloidal length element and $B_p$ the amplitude of the poloidal field. By deriving $N(\psi)$ with respect to $\psi$, we see that the local quantity that ought to be preserved in the compression is $n(dV/d\psi)$. Similar consideration regarding the conservation of entropy with a ratio of specific heat $\gamma = 5/3$ leads to the local preservation of the quantity $p(dV/d\psi)^{5/3}$.

Note that relaxing the constraint of perfect adiabaticity is of course possible and leads to transport equations such as (8) and (9) of ref. [11] (see also section 5.2 in ref. [10]). However, we will keep with perfect adiabaticity in this work.

### 2.2. Isotropy

It is natural to wonder whether the perpendicular and parallel pressure are compressed identically, or if an anisotropy is generated during the compression. This would significantly complicate the analysis. Fortunately, it was demonstrated in [16] that the equivalent of $d/dt (p\rho^{-\gamma}) = 0$ when the parallel and perpendicular pressures are considered separately is

$$ \frac{d}{dt} \left( \frac{p_{\perp}}{nB} \right) = 0 \qquad (2) $$

$$ \frac{d}{dt} \left( \frac{p_{\parallel} B^2}{n^3} \right) = 0 \qquad (3) $$

One checks easily that these relations are compatible with the scaling (1) if $p_{\perp}$ and $p_{\parallel}$ both scale as $C^5$. Therefore, no anisotropy is generated in the compression when only the lowest order of the scaling, equation (1), is considered. However, as the scaling is actually modified due to the effects derived in this paper, this conclusion must be changed. The pressure anisotropy induced by the modification of the scaling is analyzed in section 5.2.
## 2.3. Jacobian and Shafranov shift

The quantity $d\mathcal{V}/d\psi$ is simply the jacobian integrated over a flux surface. When we are interested in the magnetic axis, where the plasma is hottest, it is simply the jacobian itself (up to a $4\pi^2$ factor). It is obvious that the jacobian should play an important role, since it represents the volume element. If the volume element between two adjacent flux surfaces increases relatively to the FY scaling, then the density shall decrease relatively to that scaling, in order to conserve the total number of particles. This brings us to a first obvious effect of the Shafranov shift. During the compression, in first approximation, $\beta$ increases linearly with $C$ because the kinetic pressure increases faster than the magnetic pressure. As a result, the Shafranov shift of the magnetic axis increases, which has the effect of pushing it toward larger major radii. Since the jacobian of a torus contains a factor proportional to the major radius, this has the effect of decreasing the impact of the compression at the axis. Relatively to the FY scaling, the volume element increases, so that the density increases less than predicted by the scaling, which does not consider the internal adjustment of the equilibrium. The effect on pressure is similar, but amplified by a factor of 5/3.

There is a second, less obvious effect, that will be demonstrated in the following section. Assuming perfectly concentric flux surfaces and parameterizing the flux surface with their minor radius $r$ and geometrical angle $\theta$, the Jacobian is $J = Rr(d\psi/dr)^{-1}$. We see that in addition to the contribution of the magnetic axis expanding to larger major radii, there might be, during the compression, a relative modification of the $\psi(r)$ relation, so that $r(d\psi/dr)^{-1}$ is also changed relative to the FY scaling. The interpretation of this effect can be understood as follows. We expect $\psi \sim r^2$, hence $r(d\psi/dr)^{-1}$ independent of $r$. However If $\psi$ increases less rapidly with $r$, this means that a given $\psi = C^{\text{te}}$ surface has a larger volume, reducing the increase in density. Again, the effect on pressure is the same, but amplified by a factor of 5/3. In the case of this second effect, it is not obvious before doing the analysis that the effect is related to the Shafranov shift. However, we will prove that the dominant effect can indeed also be interpreted in terms of the Shafranov shift, and that it doubles the contribution of the first effect.

The analytical results are presented in details in the next section.

# 3. Analytical results

## 3.1. Geometry, normalization and ansatz for the fields

Before tackling the problem, we normalize the quantities so that we extract all the spatial dependence of the scaling equilibrium. We assume we start from a large radius, low magnetic field and negligible (but not vanishing) pressure equilibrium, called “standard” equilibrium. We note $B_0^\circ$ the magnetic field amplitude at the magnetic axis for the standard equilibrium, and $R_0^\circ = (\max(R) + \min(R))/2$, that is, $R_0^\circ$ is the center of the circle defining the boundary, and may differ from the major radius at the magnetic axis even for the standard equilibrium, because of a possible non-vanishing Shafranov shift. During the compression, the external boundary is reduced homothetically, and each equilibrium is parameterized with $R_0 = R_0^\circ/C$,
where $C > 1$ is the linear compression factor. The inverse aspect ratio $\epsilon = a/R_0 \ll 1$, where $a$ is the minor radius of the boundary, is kept constant. In principle, even for a circular boundary, at high pressure we should consider an elliptical deformation of the inner surfaces [5]. However, the effect is small and found not to affect our results, so it will be neglected. Thus, we assume the inner surfaces to remain perfectly circular, with a well defined radius $r$. We define a normalized minor radius $\tilde{r}$ such that $r = \epsilon \tilde{r} R_0^\circ / C$. This ensures that for all equilibria, $\tilde{r}$ varies from 0 at the axis to 1 at the boundary. Allowing for the Shafranov shift, ordered as $\Delta(r) = \epsilon^2 R_0 \tilde{\Delta}(\tilde{r})$, we can write the major radius $R$ and vertical position $Z$ in terms of $\tilde{r}$ and a geometrical angle $\theta$ as

$$R = \frac{R_0^\circ}{C} (1 + \epsilon \tilde{r} \cos \theta + \epsilon^2 \tilde{\Delta}(\tilde{r})) \quad (4)$$

$$Z = \frac{R_0^\circ}{C} \epsilon \tilde{r} \sin \theta \quad (5)$$

Note that our convention for the Shafranov shift is the same as in the Wesson book, Tokamaks [17]. The displacement of the axis is $\Delta(0)$, while the displacement of the last flux surface in $r = a$ vanishes, and we also have $\Delta'(0) = 0$. As explained in appendix C, the use of a different convention where $\Delta(0) = 0$ does not prevent from doing calculations, but prevents one from reaching the relevant conclusion. That is because the position of the magnetic axis then takes the role of an input parameter, while it is actually the result of an adjustment of the internal profiles to the external boundary conditions imposed by the wall.

To avoid any ambiguity, we will not remove the tildes on the normalized quantities. The absence of tilde indicates a quantity in MKS units, or without dimension. The jacobian is then given by

$$\begin{aligned} \mathcal{J} &= \frac{R}{d\psi/dr} (\partial_r R \partial_\theta Z - \partial_\theta R \partial_r Z) \\ &= \frac{\epsilon^2 R_0^\circ}{C^3 B_0^\circ} \frac{\tilde{r}}{\tilde{\psi}'} (1 + \epsilon \tilde{r} \cos \theta + \epsilon^2 \tilde{\Delta}) (1 + \epsilon \tilde{\Delta}') \end{aligned} \quad (6)$$

where a prime denotes derivation with respect to $\tilde{r}$. The factor $1 + \epsilon \tilde{\Delta}'$ comes from the non-orthogonality of $\nabla r$ and $\nabla \theta$ when the Shafranov shift depends on $r$ (if it did not depend on $r$, the Shafranov shift would only be a relabelling of the surfaces).

In the following, we will carry out a small parameter expansion of Grad-Shafranov's equation,

$$\Delta^* \psi = -\mu_0 R^2 \frac{dp}{d\psi} - T \frac{dT}{d\psi} \quad (7)$$

where $\Delta^* \psi \equiv R^2 \nabla \cdot (R^{-2} \nabla \psi)$ and $T = R^2 \mathbf{B} \cdot \nabla \varphi$. In the small parameter expansion, the small parameter is related to $\beta$. The idea is that when the compression starts, $\beta$ is very small, so that the influence of the pressure term in the Grad Shafranov equation is negligible compared to the current term $-TdT/d\psi$. As the compression goes on, $\beta$ increases roughly linearly with $C$ (our goal is actually to study the correction with respect to that scaling), so that the influence of the pressure in the Grad Shafranov equation increases. However, $\beta$ is not an externally fixed parameter. Externally, the operator can only fix the position of the liquid wall
surface, and the profiles will adjust in accordance with the Grad Shafranov equation. In our framework where the aspect ratio and shape of the boundary is fixed, this corresponds to using $C$ as the externally imposed parameter. As $C$ increases, the influence of the pressure becomes stronger in the Grad Shafranov equation. Therefore, it is natural to expand the pressure as $p = p_0 C^5 + p_1 C^6$, that is, the pressure scaling has a correction term that becomes increasingly significant when $C$ increases. Note that $C = 1$ corresponds to our standard equilibrium (with negligible but non-zero pressure), while $C = 0$ corresponds to zero pressure, which comes with infinite radius and vanishing magnetic field; the compressed case has $C \gg 1$. In the numerical case used to compare with the analytical results, the maximum compression has $C = 10000$. This is only a matter of convention. Since the $C = 0$ case cannot be reached in the simulations, we simply go to large radius for a given thermal energy content, which effectively decreases the pressure influence. We set the large radius case as $C = 1$ by convention. The deviation to the scaling starts to be important when $\beta$ reaches a few percents. This happens at a value of $C$ which can be large if our initial state has a very large radius. This is only because we want to compare the analytical results with the code in the case of a circular boundary. In an actual experiment, the compression ratio would not exceed about 10 [2]. Since, at lowest order, $\beta \propto C$, the maximum value reached by $C$ reflects the ratio between the final and initial $\beta$. The fact that $C = 1$ has negligible pressure means that for $C = 1$, we have $p^\circ = p_0 + p_1 \simeq p_0$, which means that we can confuse $C = 0$ and $C = 1$. We will use this several times in the following analysis. Finally, also note that mathematically, there is no requirement on the expansion parameter of a series to be small. The only thing that matters is the convergence radius of a series, and there is no reason for this radius to be smaller than 1. We now give the ordering of all the quantities used in the derivation.

The poloidal flux is normalized with $R_0^2 B_0^\circ$, so we write it as

$$ \psi(r) = R_0^2 B_0^\circ (\tilde{\psi}_0(\tilde{r}) + C \tilde{\psi}_1(\tilde{r}) + \mathcal{O}(C^2)). \quad (8) $$

Note the presence of $r$ in meters in the left hand side, and of the normalized radius $\tilde{r}$ in the right hand side. Our convention is to evaluate dimensioned quantities on dimensioned variables, and normalized quantities on normalized variables. Now, it is important to realize that the invariance of $\psi$ does not imply the invariance of $\psi(r)$. Indeed, the flux surfaces can be deformed. However, the conservation of flux does impose boundary conditions to $\tilde{\psi}_1$, namely that $\tilde{\psi}_1(0) = \tilde{\psi}_1(1) = 0$. Here, our convention that $r$ is the radius of a flux surface, and not the distance to a specific point defining the coordinate system, is essential.

The toroidal magnetic function, $T = R^2 \mathbf{B} \cdot \nabla \varphi$ is roughly proportional to $C$, since $R \propto C^{-1}$ while $B_t \propto C^2$, with $B_t$ the toroidal field amplitude. Thus, in the limit of $C \to 0$, $T$ tends to zero, but $T dT/d\psi$ remains commensurable with $\Delta^* \psi$ (while the pressure term becomes negligible). For $T$, we take

$$ T(r) = R_0^2 B_0^\circ C (\tilde{T}_0(\tilde{r}) + C \tilde{T}_1(\tilde{r}) + \mathcal{O}(C^2)). \quad (9) $$

The basic pressure scaling comes from the local invariance of $p(d\mathcal{V}/d\psi)^{5/3}$. Since $(d\mathcal{V}/d\psi) \propto C^{-3}$, we have that $p \propto C^5$. So it is natural to normalize and scale the pressure as

$$ p(r) = \frac{B_0^\circ C^5}{\mu_0} (\tilde{p}_0(\tilde{r}) + C \tilde{p}_1(\tilde{r}) + \mathcal{O}(C^2)). \quad (10) $$
The term $\tilde{p}_1$ will not enter in the calculation but be computed a posteriori, in terms of $\tilde{\psi}_1$ and other perturbed quantities. This will be our final result, allowing to estimate the modification of axial pressure with respect to FY scaling, measured by the amplitude of $\tilde{p}_1(0)$.

Finally, the invariance of $q(\psi)$ while $\psi(r)$ changes means that $q(r)$ is also not invariant. It is given by

$$q(r) = q_0(\tilde{r}) + C \frac{q'_0}{\tilde{\psi}'_0} \tilde{\psi}_1(\tilde{r}) + \mathcal{O}(C^2). \qquad (11)$$

In the following, we will use the well-known relation between $T$ and $q$:

$$T = \frac{2\pi q}{\int d\theta \mathcal{J} R^{-2}} \qquad (12)$$

The Shafranov shift is also expanded in terms of $C$, since there is a part related to the current, and a part related to the pressure, that will be negligible when $C \to 0$:

$$\tilde{\Delta} = \tilde{\Delta}_0 + C\tilde{\Delta}_1 + \mathcal{O}(C^2). \qquad (13)$$

Pressure, poloidal flux and $T$ are related through the Grad-Shafranov equation. In our framework, we are only interested in the radial dependence of $\psi$, and note that $r$ is not the radial distance to the magnetic axis (which would change on a flux surface), but the radius of the flux surfaces. So we only need to consider the lowest order Grad-Shafranov equation in terms of an expansion in powers of $\epsilon$, combined with the Shafranov shift equation [10, 17], that is:

$$\frac{1}{r} \frac{d}{dr} \left( r \frac{d\psi}{dr} \right) = -\mu_0 R_0^2 \frac{dp}{d\psi} - T \frac{dT}{d\psi} \qquad (14)$$

$$\frac{d}{dr} \left( r \left( \frac{d\psi}{dr} \right)^2 \frac{d\Delta}{dr} \right) = \frac{r}{R_0} \left( 2\mu_0 r \frac{dp}{dr} - \frac{1}{R_0^2} \left( \frac{d\psi}{dr} \right)^2 \right). \qquad (15)$$

## 3.2. Pressure deviation with respect to FY scaling

We can now expand equations (12), (14) and (15) by plugging in equations (4), (8), (9), (10), (11), and (13) and equalizing order by order in terms of the expansion parameter $C$. Our normalization allows to cancel out all possible $B_0^\circ$ and $R_0^\circ$ factors. Transforming flux derivatives $d/d\psi$ into $(d\psi/dr)^{-1}d/dr$, the Grad-Shafranov equation expands, at the two lowest orders, into:

$$\tilde{\psi}'_0 \left( \tilde{\psi}''_0 + \frac{\tilde{\psi}'_0}{\tilde{r}} \right) = -\epsilon^2 \tilde{T}_0 \tilde{T}'_0 \qquad (16)$$

$$\tilde{\psi}'_0 \left( \tilde{\psi}''_1 + \frac{\tilde{\psi}'_1}{\tilde{r}} \right) + \tilde{\psi}'_1 \left( \tilde{\psi}''_0 + \frac{\tilde{\psi}'_0}{\tilde{r}} \right) = -\epsilon^2 \left( (\tilde{T}_0 \tilde{T}_1)' + \tilde{p}'_0 \right) \qquad (17)$$

As expected, the lowest order does not depend on pressure, and the next order gives $\tilde{\psi}_1$ as a function of the unperturbed pressure profile $\tilde{p}_0$ and $\tilde{T}_1$. $\tilde{p}_0$ can be considered to be given, since it represents the shape of the pressure profile at the beginning of the compression. However, we need a relation between $\tilde{\psi}_1$ and $\tilde{T}_1$. It is provided by equation (12). We first need to analyze $\int d\theta \mathcal{J} R^{-2}$. It is given by (see Appendix A)

$$\left( \int \mathcal{J} R^2 d\theta \right)^{-1} = \frac{1}{\tilde{r}} \tilde{\psi}' \frac{CB_0^\circ R_0^\circ}{2\pi\epsilon^2} \left[ 1 - \epsilon^2 \left( \frac{\tilde{r}^2}{2} - \tilde{\Delta} - \frac{1}{2}\tilde{r}\tilde{\Delta}' \right) + \mathcal{O}(\epsilon^4) \right]. \qquad (18)$$
We rename the last square bracketed factor $f$. Owing to the expansion of $\tilde{\Delta}$ in terms of $C$, we also have $f = f_0 + C f_1 + \mathcal{O}(C^2)$, with

$$f_0 = 1 - \epsilon^2 \left( \frac{\tilde{r}^2}{2} - \tilde{\Delta}_0 - \frac{1}{2} \tilde{r} \tilde{\Delta}'_0 \right) \qquad (19)$$

$$f_1 = \epsilon^2 \left( \tilde{\Delta}_1 + \frac{1}{2} \tilde{r} \tilde{\Delta}'_1 \right). \qquad (20)$$

Replacing (18) in (12), we obtain $(\tilde{T}_0 + C\tilde{T}_1) = \frac{1}{\epsilon^2 \tilde{r}} (\tilde{\psi}'_0 + C\tilde{\psi}'_1) (f_0 + C f_1) (q_0 + C q_1)$, that is, at the two lowest orders:

$$\tilde{T}_0 = f_0 \frac{q_0}{\epsilon^2 \tilde{r}} \frac{d\tilde{\psi}_0}{d\tilde{r}} \qquad (21)$$

$$\tilde{T}_1 = f_0 \frac{1}{\epsilon^2 \tilde{r}} \frac{d}{d\tilde{r}} [q_0 \tilde{\psi}_1] + \tilde{T}_0 \frac{f_1}{f_0}, \qquad (22)$$

where $q_1 = (q'_0/\tilde{\psi}'_0)\tilde{\psi}_1$ has been used. We now replace this expression for $\tilde{T}_1$ in (17), also using (16) and (21) to eliminate $\tilde{\psi}'_0$ and $\tilde{\psi}''_0$. After multiplying by $\tilde{r}^3 (q_0 f_0 \tilde{T}_0)^{-1}$, we obtain

$$\tilde{r}^2 \tilde{\psi}''_1 \left[ 1 + \epsilon^2 \frac{\tilde{r}^2}{q_0^2 f_0^2} \right] - \tilde{r} \tilde{\psi}'_1 \left[ 1 - 2\bar{s} - \epsilon^2 \frac{\tilde{r}^2}{q_0^2 f_0^2} - \tilde{r} \frac{f_0'}{f_0} \right] + \tilde{\psi}_1 \left[ \bar{s} \left( \tilde{r} \frac{\tilde{T}_0'}{\tilde{T}_0} + \tilde{r} \frac{f_0'}{f_0} - 1 \right) + \tilde{r}^2 \frac{q_0''}{q_0} \right] = - \frac{\epsilon^2 \tilde{r}^3}{q_0 f_0 \tilde{T}_0} \frac{d}{d\tilde{r}} \left( \tilde{p}_0 + \frac{\tilde{T}_0^2}{f_0} f_1 \right). \qquad (23)$$

with $\bar{s} \equiv \tilde{r}q'/q$ the magnetic shear. Since both $\tilde{T}_0$ and $f_0$ are equal to one plus a term of order $\epsilon^2$, the above differential equation can be considerably simplified by neglecting all $\epsilon^2$ corrections. This yields the model equation

$$\tilde{r}^2 \tilde{\psi}''_1 - \tilde{r} \tilde{\psi}'_1 (1 - 2\bar{s}) - \tilde{\psi}_1 \left[ \bar{s} - \frac{\tilde{r}^2 q_0''}{q_0} \right] = - \frac{\epsilon^2 \tilde{r}^3}{q_0} \frac{d}{d\tilde{r}} (\tilde{p}_0 + f_1). \qquad (24)$$

This equation, solved with the boundary conditions $\tilde{\psi}_0(0) = \tilde{\psi}_0(1) = 0$, yields the sought-after correction to the flux-radius relation. We only need to express the $f_1$ term, which comes from the Shafranov shift equation (15). At the two lowest orders, it yields

$$\frac{d}{d\tilde{r}} \left[ \tilde{r} \left( \tilde{\psi}'_0 \right)^2 \frac{d\tilde{\Delta}_0}{d\tilde{r}} \right] = -\tilde{r} \left( \tilde{\psi}'_0 \right)^2 \qquad (25)$$

$$\frac{d}{d\tilde{r}} \left[ \tilde{r} \left( \tilde{\psi}'_0 \right)^2 \frac{d\tilde{\Delta}_1}{d\tilde{r}} \right] = 2\epsilon^2 \tilde{r}^2 \tilde{p}'_0 - 2 \frac{d}{d\tilde{r}} \left[ \tilde{r} \tilde{\psi}'_0 \tilde{\psi}'_1 \frac{d\tilde{\Delta}_0}{d\tilde{r}} \right] - 2\tilde{r} \tilde{\psi}'_0 \tilde{\psi}'_1. \qquad (26)$$

It can be proven that the terms containing $\tilde{\psi}'_1$ lead only to $\epsilon^2$ contributions (see Appendix B), which we neglect. By replacing $\tilde{\psi}'_0$ using expression (21) and replacing $\tilde{T}_0^2/f_0^2$ with 1 yields the model equation for $\tilde{\Delta}_1$:

$$\frac{d}{d\tilde{r}} \left( \epsilon^2 \frac{\tilde{r}^3}{q_0^2} \frac{d\tilde{\Delta}_1}{d\tilde{r}} \right) = 2\tilde{r}^2 \tilde{p}'_0. \qquad (27)$$

Since it is a Shafranov shift equation, it ought to be solved with boundary conditions $\tilde{\Delta}_1(1) = 0$ and $\tilde{\Delta}'_1(0) = 0$. Since $f_1 \sim \epsilon^2 \tilde{\Delta}_1$, equation (27) shows that both $\tilde{p}_0$ and $f_1$ terms in equation (24) are of the same order of magnitude.

This clarifies the interpretation of the phenomenon. When the compression ramps up, the pressure term kicks in the Grad Shafranov equation (14), yielding a contribution to the
modification of the flux-radius relation. However, it also leads to an increase of the Shafranov shift, which gives an additional contribution to the Grad Shafranov equation through the $f_1$ term. In fact, we find that the latter dominates the former contribution (the $\tilde{p}_0$ term in the right hand side of equation (24)) by a factor of roughly $q^2$ (see section 3.3).

In order to derive the pressure perturbation, there is one last preliminary quantity, which is $d\mathcal{V}/d\psi$. Poloidal integration of (6) yields

$$ \frac{d\mathcal{V}}{d\psi} = \frac{4\pi^2\varepsilon^2 R_0^\circ \tilde{r}}{C^3 B_0^\circ \tilde{\psi}'} \left[ 1 + \varepsilon^2 \left( \tilde{\Delta} + \frac{1}{2}\tilde{r}\tilde{\Delta}' \right) \right] \quad (28) $$

The bracketed factor can be written $e_0 + C f_1 + \mathcal{O}(C^2)$, with $e_0 = 1 + \varepsilon^2 (\tilde{\Delta}_0 + \frac{1}{2}\tilde{r}\tilde{\Delta}'_0) = f_0 + \varepsilon^2 \tilde{r}^2 / 2$.

The invariance of $p(d\mathcal{V}/d\psi)^{5/3}$ means that the pressure at a given compression ratio $C$ is given by

$$ p(\psi) = p^{\circ}(\psi) \left( \frac{d\mathcal{V}^{\circ}/d\psi}{d\mathcal{V}/d\psi} \right)^{5/3}, \quad (29) $$

where $\circ$ stands for the standard equilibrium, where $C=1$. Recalling the standard equilibrium is defined as one for which the effect of the pressure is negligible, we must have all quantities with index 1 much smaller than one. With this approximation (which amounts to confuse $C=0$ and $C=1$), we can expand equation (29) into

$$ (\tilde{p}_0(\tilde{r}) + C \tilde{p}_1(\tilde{r})) = \left( \tilde{p}_0 + C \frac{\tilde{p}'_0}{\tilde{\psi}'_0} \tilde{\psi}_1 \right) (1 + C f_1)^{-5/3} \left( 1 + C \frac{\tilde{\psi}'_1}{\tilde{\psi}'_0} \right)^{5/3}, \quad (30) $$

which yields the perturbation of the pressure profile with respect to the FY scaling, $\delta p/p = C\tilde{p}_1/\tilde{p}_0$, measured by the following quantity:

$$ \frac{\tilde{p}_1}{\tilde{p}_0} = -\frac{5}{3} \left( f_1 - \frac{\tilde{\psi}'_1}{\tilde{\psi}'_0} \right) + \frac{\tilde{p}'_0}{\tilde{\psi}'_0} \frac{\tilde{\psi}_1}{\tilde{p}_0}. \quad (31) $$

At the axis, since $\tilde{\psi}_1(0) = 0$, the last term vanishes, and we see the pressure deviation is indeed the sum of two terms. At the axis, $f_1 = \varepsilon^2 \tilde{\Delta}_1$, and the $-5/3 f_1$ term accounts for the reduction due to the magnetic axis being displaced to larger major radius regions with larger volume element, the first effect described in section 2.3. Note that $\tilde{\Delta}_1$ is the contribution to the Shafranov shift of the pressure. Since $\beta$ increases with $C$, so does $\tilde{\Delta}$, so that $\tilde{\Delta}_1$ is a positive quantity. The second effect is accounted for by the $5/3 \tilde{\psi}'_1/\tilde{\psi}'_0$ term. As shown in section 3.3, this term is almost equal to and adds to the first one. It is also an indirect effect of the shift because $\tilde{\psi}'_1$ is mainly driven by $f_1$ in equation (24).

As noted in the introduction, Greene, Johnson and Weimer have already derived a correction to the Furth-Yoshikawa scaling. [5]. However, although the result we find bears similarities with theirs, their expression cannot be directly exploited. The main reason lies in a convention chosen for the compression parameter, which amounts to assuming that the position of the magnetic axis itself is controlled, while only the external boundary conditions are experimentally controlled. More details are presented in Appendix C.
### 3.3. The case of vanishing shear

We can gain a lot of insight by computing exactly the effect in the case of vanishing shear. The existence of a sheared safety factor profile will affect the numbers and prevent the derivation of a closed form, but does not change the nature of the effect. In this section, we simply assume that the safety factor is constant and we denote it simply by $q$. Then, we are solving

$$ \tilde{r}^2 \tilde{\psi}_1'' - \tilde{r} \tilde{\psi}_1' = - \frac{\varepsilon^2 \tilde{r}^3}{q} (\tilde{p}_0 + f_1)' \quad (32) $$

$$ (\varepsilon^2 \tilde{r}^3 \tilde{\Delta}_1)' = 2 \tilde{r}^2 q^2 \tilde{p}_0' \quad (33) $$

$$ f_1 = \varepsilon^2 \left( \tilde{\Delta}_1 + \frac{1}{2} \tilde{r} \tilde{\Delta}_1' \right). \quad (34) $$

We shall be interested in the effect at the magnetic axis only, because it can be related to global quantities, and because it is where the plasma is hottest. We introduce $\beta^\circ \equiv 2\mu_0 p^\circ / B_0^{\circ 2}$ and $\beta_P^\circ = 2\mu_0 \overline{p^\circ} / B_\theta^{\circ 2}$, respectively the axis and global poloidal beta of the standard equilibrium, with $B_\theta^\circ$ taken at the last flux surface and $\overline{p^\circ}$ the volume averaged pressure. In our normalization, using $B_\theta^\circ(a^\circ)/B_0^\circ = \varepsilon/q$, valid for small aspect ratio, this is:

$$ \beta^\circ = 2 \tilde{p}_0(0) \quad (35) $$

$$ \beta_P^\circ = \frac{4q^2}{\varepsilon^2} \int_0^1 \tilde{p}_0(\tilde{r}) \tilde{r} d\tilde{r} \quad (36) $$

First, we compute the $f_1$ contribution at the axis, which is simply $\varepsilon^2 \tilde{\Delta}_1(0)$. This yields (see Appendix D):

$$ f_1(0) = \frac{q^2}{2} \left( \beta^\circ - \frac{\varepsilon^2}{q^2} \beta_P^\circ \right). \quad (37) $$

Note that $\varepsilon^2 q^{-2} \beta_P^\circ$ is smaller than $\beta^\circ$, making the Shafranov shift a positive quantity, as is well known. For instance for a parabolic pressure profile, $f_1(0) = \varepsilon^2 \tilde{\Delta}_1 = q^2 \beta^\circ / 4$.

To compute the contribution of $\tilde{\psi}_1'/\tilde{\psi}_0'$, we note that since typically, the poloidal flux varies as $\tilde{r}^2$ close to the axis, we need to look at the ratio $(\tilde{\psi}_1'/\tilde{r})/(\tilde{\psi}_0'/\tilde{r})$. Equation (32) is an Euler-Cauchy equation with solution

$$ \tilde{\psi}_1(\tilde{r}) = K \frac{\tilde{r}^2}{2} - \frac{\varepsilon^2}{q} \int_0^{\tilde{r}} ds \, s(\tilde{p}(0) - \tilde{p}(s) + f_1(0) - f_1(s)), \quad (38) $$

where $K = (\tilde{\psi}_1'/\tilde{r})|_{\tilde{r}=0}$ is an integration constant such that $\tilde{\psi}_1(0) = 0$. Therefore, eliminating $\tilde{\psi}_0'$ using equation (21) with $\tilde{T}_0/f_0$ replaced with 1 (valid up to order $\varepsilon^2$ terms), we find

$$ \begin{aligned} \frac{\tilde{\psi}_1'/\tilde{r}}{\tilde{\psi}_0'/\tilde{r}} &= \frac{qK}{\varepsilon^2} \\ &= -2 \int_0^1 d\tilde{r} \, \tilde{r} (\tilde{p}(0) - \tilde{p}(\tilde{r}) + f_1(0) - f_1(\tilde{r})) \\ &= -\frac{1}{2} \left( \beta^\circ - \frac{\varepsilon^2}{q^2} \beta_P^\circ \right) (1+q^2), \end{aligned} \quad (39) $$
where we have used equation (37) as well as the fact that $\int_0^1 \tilde{r} f_1(\tilde{r}) d\tilde{r} = 0$ by integration by parts of (34).

We can combine (39) and (37) into (31) at the axis to find the pressure perturbation. By noting that at lowest order in $C$, the FY scaling applies, $\beta \simeq C\beta^\circ$, $\beta_P \simeq C\beta_P^\circ$, we finally find:

$$ \left. \frac{\delta p}{p} \right|_{r=0} = -\frac{5}{6} \left( \beta - \frac{\varepsilon^2}{q^2} \beta_P \right) (1 + 2q^2). \qquad (40) $$

The $\delta$ here is to be interpreted as “with respect to the FY scaling”. The first important conclusion that can be drawn from this result is that the deviation is always negative, meanign that the pressure increases less than predicted by the FY scaling. Furthermore, the analytical result can be seen as being made of two contributions, with relative amplitudes respectively 1 and $2q^2$, where the latter can be clearly decomposed into two $q^2$ contributions. Since in a tokamak a large part of the plasma has $q > 1$, the $2q^2$ contribution is dominant. The first contribution (with amplitude 1) comes from the $\tilde{p}_0$ term on the right hand side of (24). The remaining two $q^2$ contributions come respectively from the two effects advertized in section 2.3: (i) the Shafranov shift takes the axis to larger major radii, and (ii) it also modifies the flux vs radius relation by opening the flux surfaces, through the $f_1$ term in equation (24). As one can immediately see, the effect is not small if $\beta$ reaches several tens of percents, as is planned in General Fusion concepts [2].

# 4. Comparison with CHEASE simulations

## 4.1. Setting a predefined safety factor profile in CHEASE and enforcing entropy conservation

The CHEASE code [15] solves Grad Shafranov equation,

$$ R^2 \nabla \cdot \left( \frac{1}{R^2} \nabla \psi \right) = -\mu_0 R^2 \frac{dp}{d\psi} - T \frac{dT}{d\psi} \qquad (41) $$

with a prescribed boundary, and several options to implement the $p(\psi)$ and $T(\psi)$ profiles. Although $p'$ is reasonably constrained from the experimental knowledge of the radial pressure profile, $T$ is more difficult to get right, and arbitrary choices will often lead to unrealistic current or safety factor profiles. It is, therefore, more usual to choose the current profile, and to use the relation between this current and $T$ to define the profile of the latter. Namely, a common choice is to define the surface averaged current density as

$$ I^*(\psi) = \frac{\oint (\mathbf{j} \cdot \nabla \varphi) \mathcal{J} d\theta}{\oint R^{-1} \mathcal{J} d\theta} = \frac{\langle j_\varphi / R \rangle}{\langle 1/R \rangle}, \qquad (42) $$

where the flux surface average is defined by

$$ \langle A \rangle = \left( \oint \mathcal{J} d\theta \right)^{-1} \oint A \mathcal{J} d\theta. \qquad (43) $$

This current is related to $dp/d\psi$ and $TdT/d\psi$ via the following relation [15]:

$$ I^* = - \frac{C_1}{C_0} \frac{dp}{d\psi} - \frac{C_2}{C_0} T \frac{dT}{d\psi}, \qquad (44) $$
where

$$ \{C_0, C_1, C_2\} = \oint \left\{ \frac{1}{R}, 1, \frac{1}{R^2} \right\} \mathcal{J} d\theta. \qquad (45) $$

Since the integrals (45), carried at $\psi = C^{\text{te}}$, obviously depend on the $\psi$ solution, the relation (44) is iterated over, every time a new solution is computed, until convergence is achieved.

In some cases, and for obvious reasons in the present case, it is desirable to define $T$ not through the current, but directly by inputting the safety factor profile $q(\psi)$. The procedure, which is non trivial, was implemented in CHEASE in 2016 and is described in reference [18]. The procedure is now outlined.

In principle, fixing the safety profile could be done by using equation (12), which is simply $T = 2\pi q/C_2$ with the newly introduced notation. Unfortunately, this does not work as well as for the current. If one attempts to input $I^*$ in the form

$$ I^* = -\frac{C_1}{C_0} \frac{dp}{d\psi} - \frac{4\pi^2}{C_0} q \frac{d}{d\psi} \left( \frac{q}{C_2} \right), \qquad (46) $$

the safety profile rapidly develops spurious oscillations instead of converging to the desired value. If one allows for a gradual iterative procedure

$$ I_{\text{input},n}^* = \alpha I_{\text{target}}^* + (1 - \alpha) I_{n-1}^* \qquad (47) $$

with $I_k^*$ the current at step $k$ and $I_{\text{target}}^*$ given by equation (46), the oscillations do not disappear even for $\alpha \ll 1$. A similar problem was already mentioned in past works [19, 11] where it was also found that direct use of equation (12) was not adequate. Instead, alternative relations between $q$, $dp/d\psi$ and $TdT/d\psi$ had to be put to use.

The solution developed for CHEASE relies on the same idea of using a more sophisticated relation between these quantities. The total current $I_p$ contained within a constant $\psi$ surface can be obtained either as a surface integral of the current, or through Ampère's law:

$$ \begin{aligned} I_p(\psi) &= \int^\psi d\psi \oint (\mathbf{j} \cdot \nabla \varphi) \mathcal{J} d\theta \\ &= \mu_0^{-1} \oint B_P dl \\ &= (2\pi\mu_0)^{-1} \frac{d\mathcal{V}}{d\psi} \langle B_P^2 \rangle. \end{aligned} \qquad (48) $$

In (48), $B_P^{-1} dl = \mathcal{J} d\theta$ was used. We see by taking the $\psi$ derivative of $I_p$ that an alternative definition of $I^*$ is

$$ I^* = \frac{1}{\mu_0 \langle 1/R \rangle} \frac{d}{d\mathcal{V}} \left( \langle B_P^2 \rangle \frac{d\mathcal{V}}{d\psi} \right). \qquad (49) $$

We can manipulate this expression by defining an effective radius $\rho_T$ such that the toroidal flux is $\psi_T = \rho_T^2/2$, with $q = d\psi_T/d\psi$. This means that $d\psi/d\rho_T = q^{-1}\rho_T$. Then, by also using $B_P = |\nabla\psi|/R$, we can reexpress the term in the brackets of $I^*$ to obtain

$$ I^* = \frac{1}{\mu_0 \langle 1/R \rangle} \frac{d}{d\mathcal{V}} \left( \frac{d\mathcal{V}}{d\rho_T} \left\langle \frac{|\nabla \rho_T|^2}{R^2} \right\rangle \frac{\rho_T}{q} \right). \qquad (50) $$
The apparition of $q$ at the denominator suggests the following procedure. Denoting by $q_{\text{target}}$ the safety factor one wishes to obtain in the final equilibrium, and $q$ the safety factor corresponding to the currently available $\psi$ solution (or first guess), one replaces equation (49) with

$$I^* = \frac{1}{\mu_0 \langle 1/R \rangle} \frac{d}{d\mathcal{V}} \left( \langle B_P^2 \rangle \frac{d\mathcal{V}}{d\psi} \frac{q}{q_{\text{target}}} \right). \quad (51)$$

This leads to the disappearance of the oscillations and to very good convergence properties. Note that equation (47) is also used. A value of $\alpha = 0.2$ has given good results in our case. The iterative procedure of inputting $I^*$ through (51) leads to an equilibrium with the desired safety factor. More details regarding the implementation can be found in reference [18].

The question that remains is how to ensure the invariance of $p(d\mathcal{V}/d\psi)^{5/3}$. Fortunately, there is no difficulty here. Starting from an initial equilibrium (say, our standard equilibrium), we first update the pressure profile according to the FY scaling to obtain an initial guess, used to compute an equilibrium. At the next iteration, the pressure profile is defined simply as

$$p = \frac{\left[ p \left( \frac{d\mathcal{V}}{d\psi} \right)^{5/3} \right]_{\text{target}}}{\left( \frac{d\mathcal{V}}{d\psi} \right)^{5/3}}, \quad (52)$$

where at the numerator the quantity is fixed once and for all, and the denominator is computed from the last equilibrium. The operation is repeated as many times as necessary. An $\mathcal{L}^2$ convergence to a tolerance of $10^{-6}$ can be obtained usually in roughly half a dozen iterations. The simplified CHEASE workflow adapted to the generation of flux and entropy conserving equilibria is represented in figure 4.1

<figure><img src="image_2.png" /><figcaption>Figure 1. CHEASE workflow for the generation of flux and entropy conserving equilibria.</figcaption></figure>

## 4.2. Results with a circular boundary

In this section, the analytical results of section 3 are tested against CHEASE computations. In order to be as close as possible to the theory, we impose a circular boundary, and use a very
<figure><img src="image_3.png" /><figcaption>Figure 2. Comparison between the analytical and numerical profiles of (a) $\tilde{\psi}_1$ and (b) $\tilde{p}_1$.</figcaption></figure>

large aspect ratio $\varepsilon^{-1} = 10$. In order to make the pressure of the standard equilibrium ($C = 1$) as negligible as possible, we choose $R_0^\circ = 1200\text{m}$ and $B_0^\circ = 7 \times 10^{-7}\text{T}$. At $C = 1000$, we obtain a plasma with the starting parameters of General Fusion [2], $R_0 = 1.2\text{m}$ and $B_0 = 0.7\text{T}$. The safety factor profile is monotonous, with $q_{r=0} = 0.8$, $q_{r=a} = 5$ and $q = 1$ roughly at mid-radius. The pressure profile is also monotonous and peaked in the center. We have generated 30 equilibria between $C = 1$ and $C = 10000$. For each equilibrium, the radius of the flux surfaces can be measured to infer the $\psi(r)$ dependence. This allows to compute $\tilde{\psi}_1$ point by point: for every flux surface with radius $\tilde{r}$, $\tilde{\psi}(\tilde{r})$ is computed as the linear coefficient of a polynomial fit of $\tilde{\psi} - \tilde{\psi}_0$, where $\tilde{\psi}_0$ is simply the normalized flux of the standard equilibrium. Note that only a small error is introduced by confusing the standard equilibrium with $C = 0$, because $|\tilde{\psi}_1| \ll |\tilde{\psi}_0|$. The numerical result is compared with the solution of equation (24), with $f_1$ computed from (27). Similarly, once $f_1$ and $\tilde{\psi}_1$ are computed, the pressure deviation $\tilde{p}_1$ can be computed and compared with its numerical counterpart, obtained with the same procedure as for $\tilde{\psi}_1$. The results can be seen in figure 2. The agreement is excellent for $\tilde{\psi}_1$ with an $L^2$ norm difference of only 1.5%, that is, of order $\varepsilon^2$ as could be expected. The agreement is also very good for $\tilde{p}_1$, with a difference lower than 5%.

Next, we show in figure 3 how well the analytical result captures the reduction of the axis pressure with respect to the FY scaling when $C$ is increased. Figures 3 a), b) and c) respectively show the axis beta, the poloidal beta and the toroidal current. We see clear deviations with respect to the FY scaling for all these quantities. Figure 3 d) shows the axial pressure divided by $C^5$. The FY scaling in this case is merely a horizontal straight line. We show in the red dashed-dotted curve that our analytical result indeed captures the deviation with respect to the scaling for modest values of $C$. For $C > 3000$, the plain and red dashed-dotted curves separate, meaning the next order in the $C$ expansion, which we have not calculated, kicks in. In principle, the next order can be computed with this theory, the algebra being quite similar. However, we have found that the Shafranov shift equation, (15), cannot
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison between the numerical result and either the FY scaling or the corrected FY scaling in the case of circular cross section, for (a) the axis beta, (b) the poloidal beta, (c) the toroidal current, (d) the axis pressure scaled down by $C^5$.</figcaption></figure>

adequately reproduce the next order of the shift, $\tilde{\Delta}_2$. This might be because the ellipticity of the surfaces starts to play an important role, indeed the derivation of the Shafranov shift equation explicitly assumes circularity of the flux surfaces.

We can, in fact, obtain an even better agreement if instead of using the Shafranov shift equation (15) to obtain the shift, we use the one measured in the CHEASE simulation. Then, $\delta p/p$ at the axis contains a first contribution $-5/3\epsilon^2(\tilde{\Delta} - \tilde{\Delta}^\circ)$, the intuitive geometric effect, and the other contribution of the shift coming through $\tilde{\psi}'_1$ at the axis. We have seen that the latter is almost exactly equal to the former in the case of zero shear, up to a term that is subdominant unless $q \simeq 1$, which is never the case over the whole volume in tokamak plasmas. This suggests that even for non-vanishing shear, it may suffice to double the geometric effect of the shift to get the whole effect. Thus, we conjecture that the following formula can be used with a good approximation:

$$ \frac{\delta p}{p} = -\frac{10}{3} \left( \frac{\Delta}{R_0} - \frac{\Delta^\circ}{R_0^\circ} \right), \qquad (53) $$

where we have removed the normalization to make it a more readily usable quantity. The disadvantage is that it assumes that the Shafranov shift is known. However, since it is an important MHD quantity, it is reasonable to assume that it may be known to operators with reasonable precision even in the absence of a complete equilibrium solution. The green dashed
curve in figure 3 d) shows that indeed, (53) provides a satisfactory approximation to describe the deviation with respect to the FY scaling over the whole C range.

4.3. Results with shaping and a tight aspect ratio

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Visualization of the flux surfaces in the case of tight aspect ratio and with shaping. The equilibria are also plotted for compression ratios $C = 2$, $C = 5$ and $C = 10$.</figcaption></figure>

The previous section served to confirm the theory by using large aspect ratio and a circular boundary. If we now add ellipticity, triangularity and use a tight aspect ratio, it is clear that the analytical theory of section 3 is out of its domain of validity. However, the geometric effect of the shift will necessarily remain the same, owing to the major radius dependence of the volume element in a torus. We can also expect that there will be an opening of the inner flux surfaces due to the compression, even if we cannot compute the effect precisely. Hence, it is tempting to test formula (53) in the case of a shaped equilibrium.

We have generated 20 equilibria with aspect ratio $\epsilon = 0.7$, elongation $\kappa = 1.3$ and triangularity $\delta = 0.4$. This looks similar to a spherical tokamak. The safety factor profile is the same as in the previous subsection. The standard equilibrium in this case ($C = 1$) is defined as having $R_0^\circ = 1.2\text{m}$ and $B_0^\circ = 0.7\text{T}$, and the maximum compression has $C = 10$. The results are seen in figure 5 and 6. Figures 5 a), b) and c) illustrate the large deviation of $\beta$ and $\beta_P$ with respect to FY, and a modest deviation for the current. Thus, the deviation of $\beta_P$ is mostly due to the pressure deviation, which is seen in 5 d). We see that, again, equation (53) captures surprisingly well the deviation, with a precision roughly similar to that observed in the case of a circular boundary.

Figure 6 a) is the same as 6 d) except that the pressure is not divided by $C^5$. This makes the deviation much more apparent. The final pressure is less than $5 \times 10^8$ Pa instead
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison between the numerical result and either the FY scaling or equation (53) in the case of D-shaped cross section with tight aspect ratio (κ = 1.3, δ = 0.4 and ε = 0.7), for (a) the axis beta, (b) the poloidal beta, (c) the toroidal current, (d) the axis pressure scaled down by C⁵.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between the numerical result and either the FY scaling or equation (53) in the case of D-shaped cross section with tight aspect ratio (κ = 1.3, δ = 0.4 and ε = 0.7) for (a) the axis pressure and (b) the fusion power.</figcaption></figure>
of almost $8 \times 10^8$ Pa predicted by the scaling. This has a strong impact on the fusion power, which, as is well known, scales roughly as $p^2$ when $T \sim 10$ keV. As seen in 6 b), where $p^2$ is represented, the effect on the fusion power is very important. The final reduction in peak fusion power is of more than 2.5, which considerably impacts the optimistic conclusions drawn in reference [2]. Note that we assume perfect adiabaticity, an assumption which will hold only if the compression time is significantly shorter than the energy and flux confinement times.

At this point, the reader may consider that this does not matter, because if the pressure increases less rapidly than anticipated, then it means that the thermal energy stored in the plasma also increases less rapidly. So the mechanical pistons performing work $-\int pdV$ will simply compress the plasma to a larger compression ratio $C$, thus cancelling the effect of the shift described in this manuscript. Unfortunately, this does not happen, and there are two reasons for that. The first reason is that in a tokamak plasma, the magnetic energy dominates, so that the thermal energy not increasing rapidly enough would have little effect on the pistons during most of the compression, except maybe at the end when $\beta$ reaches several tens of percents. The second reason is that actually, the thermal energy increases slightly more rapidly than the scaling. This is counterintuitive, but possible because, as seen in figure 2 b), the pressure perturbation is positive for larger radii, which have a much larger volume element. In any case, the resulting effect is small and the total energy computed numerically by volume integration of the CHEASE result is virtually indistinguishable from the scaling $E_{\text{tot}} = E_{\text{th}}^\circ C^2 + E_{\text{mag}}^\circ C$, where $E_{\text{th}}^\circ$ and $E_{\text{mag}}^\circ$ are respectively the thermal and magnetic energies of the standard equilibrium.

# 5. Discussion

## 5.1. Summary

We have shown that the central pressure increases much slower than predicted by the Furth-Yoshikawa scaling when the plasma is compressed adiabatically while keeping the aspect ratio constant. The reduction with respect to the scaling is of order $\beta$ and can clearly be interpreted in terms of the Shafranov shift. An effect of including the shift is to take the axis to larger major radii, where the larger volume element leads to a relative decrease of the plasma pressure. This effect is doubled by the modification of the flux vs radius relation: the inner surfaces expand relatively to the plasma boundary during the compression. Both effects are related to the Shafranov shift, as has been clearly demonstrated in section 3. The theory works surprisingly well out of its domain of validity, as seen in section 4.3. We conjecture that this would be the case also if we removed the constraint of constant aspect ratio. When the aspect ratio is not kept constant, $a$ and $R_0$ are allowed to vary independently, and the recipe to obtain the scaling is as follows. The volume scales as $V \sim a^2 R_0$, the density as $n \sim V^{-1}$, the pressure as $p \sim V^{-5/3}$, and the toroidal and poloidal fields $B_T$ and $B_P$ also scale differently because the conserved toroidal and poloidal fluxes are $\psi_T \sim a^2 B_T$ and $\psi \sim a R_0 B_P$. Nonetheless, whatever the scaling becomes, the Shafranov shift will increase if $\beta$ increases,
and this will affect the central pressure with a contribution $\delta p/p = -5/3\epsilon^2(\tilde{\Delta} - \tilde{\Delta}^\circ)$. Thus, our conjecture is that the second effect related to the flux vs radius dependence will again be present and be roughly equal to the first one, with equation (53) remaining valid. The verification is left for future work.

It this is indeed verified, then we can conclude that it is highly desirable to limit the increase in Shafranov shift in General Fusion plasmas. This can be achieved by changing the aspect ratio during the compression. The discussion is then profile dependent. Depending on the current and pressure profile, the Shafranov shift will be dominated by effects associated to $\beta$ or $\beta_p$ (as highlighted by equation (37)). When the aspect ratio is changed, $\beta$ and $\beta_p$ vary as [1]:

$$ \beta \sim \epsilon^{2/3} \frac{1}{R_0} \qquad (54) $$

$$ \beta_p \sim \epsilon^{-1/3} \frac{1}{a}, \qquad (55) $$

as can be easily verified by using the above rules. Therefore, depending on the profiles, it might be necessary to limit the Shafranov shift by trying to force the plasma to increase or decrease its aspect ratio.

Another avenue to be studied is using triangularity as an actuator. Indeed negative triangularity plasma can have a larger Shafranov shift at zero beta than positive triangularity shapes at finite beta, as discussed in ref. [20]. Therefore one could start with low or slightly negative triangularity and move towards positive triangularity at maximum compression, which would also move towards better MHD stability when reaching high beta plasmas.

## 5.2. Apparition of pressure anisotropy

In section 2.2, we have shown that the FY scaling is compatible with the hypothesis of isotropy. However, since there is a deviation with respect to that scaling, we must reevaluate this conclusion. Combining and developing equations (2) and (3), we obtain

$$ \frac{d}{dt}(p_{\perp} - p_{\|}) = \left(\frac{\dot{n}}{n} + \frac{\dot{B}}{B}\right)(p_{\perp} - p_{\|}) + \left(3\frac{\dot{B}}{B} - 2\frac{\dot{n}}{n}\right)p_{\|}, \qquad (56) $$

where the dot represents the time derivative. At lowest order, $\dot{n}/n + \dot{B}/B = 5\dot{C}/C$. So the first term in the right hand side of equation (56) tells us that the pressure anisotropy increases, like the pressure, as $C^5$. However, there has to be an anisotropy seed in the first place. The term responsible for this is the second one. At lowest order (assuming the FY scaling), we find that $3\dot{B}/B - 2\dot{n}/n = 0$, in accordance with our previous conclusion that there is no anisotropy if the FY scaling is correct. If we now consider the next order of the calculation, we can express the anisotropy generating term in terms of the Shafranov shift, as follows. To remove the $\cos\theta$ dependence of the magnetic field, let us just consider the axis. By writing $B = T/R$, we obtain, at lowest order in $\epsilon^2$ ($\tilde{T}_0 \simeq 1, 1 + \epsilon^2\tilde{\Delta}_0 \simeq 1$):

$$ B = B_0 C^2 (1 + (\tilde{T}_1 - \epsilon^2 \tilde{\Delta}_1) C). \qquad (57) $$
The density is evaluated using the formulæ for the pressure, e.g. (C.7), evaluated at the axis, and with the factor of 5/3 replaced with 1. This gives

$$n = n^{\circ} C^{3} (1 + (\tilde{T}_{1} - 2\epsilon^{2}\tilde{\Delta}_{1}) C). \qquad (58)$$

Using the latter two expressions, we find that

$$3\frac{\dot{B}}{B} - 2\frac{\dot{n}}{n} = (\tilde{T}_{1} + \epsilon^{2}\tilde{\Delta}_{1}) C\frac{\dot{C}}{C}. \qquad (59)$$

The numerical results of section 4 suggest that the $\tilde{T}_1$ term is subdominant compared to the Shafranov shift term. In that case, we can conclude that a pressure anisotropy will appear, such that $p_{\perp} > p_{\parallel}$. By integrating equation (56) for short times (and neglecting $\tilde{T}_1$), we predict an anisotropy increasing with the Shafranov shift as

$$\frac{p_{\perp} - p_{\parallel}}{p} \sim \left( \frac{\Delta}{R_0} - \frac{\Delta^{\circ}}{R_0^{\circ}} \right). \qquad (60)$$

This is small, but reaches about 10% in the case examined in section 4.3. Unfortunately, since anisotropy is not taken into account in the CHEASE code, we cannot compare this prediction with the numerical simulations. Also, note that if the anisotropy becomes too large, our analytical model breaks down because the assumption of an isotropic pressure is no longer valid.

### 5.3. Correction of other global quantities

In this paper, we have limited the discussion of the deviation with respect to the FY scaling to the pressure, and concentrated our efforts on the central pressure. In principle, we should be able to derive the first order corrections to other quantities such as $\beta$, $\beta_p$ and the toroidal currents. The deviations are observed in the plots a), b) and c) of figures 3 and 5. However, even if possible in principle, the computation of the total current, which requires radial integration of the derived expressions, does not bring much insight. Instead, for the central pressure, we have found that we can easily obtain a satisfactory ballpark estimation, with significant practical consequences. Nonetheless, we insist that the deviations of $\beta$, $\beta_p$ and $I_t$ with respect to the FY scaling do not have anything to do with the aspect ratio not being large or the plasma boundary changing shape. It is simply a result of the internal adjustment of the profiles following the Grad Shafranov equation.

### 5.4. The flux conserving tokamak

In the decade that followed the work of Furth and Yoshikawa on the scaling of tokamak adiabatic compression, the concept of flux conservation in the study of Grad-Shafranov equation remained important, albeit with a different point of view, that of the so-called flux conserving tokamak. The idea was to compute the modification of a magnetic equilibrium when an initial equilibrium is heated via neutral beam injection or wave resonant heating, which contribute only to the increase of the thermal energy of the plasma. Then, the equilibrium is modified at constant toroidal and poloidal fluxes (which implies constant safety
factor for each constant $\psi$ flux surface). At the time, when the numerical resolution of Grad-Shafranov equation was still in its infancy, the poloidal beta limitations in pure equilibrium terms were not known, for instance, it was debated whether $\beta_p$ was limited to not exceed the inverse aspect ratio $A = \varepsilon^{-1}$ [11]. The flux conserving tokamak was seen as an approach allowing high beta operation of a tokamak, and even led to reactor designs [21]. The main result of this research was the scaling of the toroidal current when the pressure is increased while conserving the magnetic flux, namely $I_t \propto p^{1/3}$ [22]. Note that this scaling cannot be directly compared to our results, because in the case of the flux conserving tokamak, the volume of the plasma does not change, or only marginally, for instance in order to reduce the skin currents that may appear [19]. Consequently, the FY paper is barely cited in the theoretical research of the flux conserving tokamak, but to underline the differences [23].

### 5.5. Connection with past numerical works

In references [11, 12, 13], numerical computations of adiabatic compression with varying aspect ratio were performed and compared with the scalings of Furth and Yoshikawa. Significant differences were found, but the reason for these differences was not well established. For instance, reference [12] attributes the difference observed during the compression toward tight aspect ratio to the fact that FY scaling is valid only for large aspect ratio. Reference [24] criticizes these conclusions by proposing another interpretation. The author considers that the compressed flux is not only the plasma flux, but also the flux trapped between the plasma boundary and the wall. His discussion then includes the surface currents arising during the compression, and attributes the difference with the scaling to these. Although the analysis is physically sound, it is not relevant to references [12, 11, 13], nor to the present manuscript, because here and in the latter references, the plasma is assumed to be in contact with the wall. Indeed, the boundary condition is that of a perfectly conducting shell, setting the shape of the external boundary of the plasma. Also, note that reference [5] is not cited in these works, despite the direct relevance of the adiabatic compression calculations performed there.

### 5.6. Surface currents

We have dismissed the explanation of reference [24] in terms of the vacuum flux between the plasma and the wall being compressed and leading to surface currents, not because it is not relevant to these plasmas, but because in our view, this does not constitute the main cause of the deviation with respect to the scaling. Indeed, we do not have to consider surface currents or an additional zone between plasma and wall to derive the result. However, we do not wish to convey the idea that the question of these surface currents is not important. One must realize the essential difference between an actual tokamak plasma and the General Fusion plasmas. The former are held in place, prevented from expanding radially under the effect of pressure and hoop forces, by an externally generated vacuum field. On the contrary, when the plasma of General Fusion expands radially, it meets the wall, where it generates image currents, eventually emulating the role of the vacuum field. This is an
essential difference of a magnetized low $\beta$ plasma compared to a hot gas in a container, which simply balances its pressure against the physical wall. When the plasma is compressed, the motion of the wall against the flux conserving plasma will generate surface currents (the inward compression motion of the plasma can be largely seen as the result of the Laplace force arising from the interaction between these currents at the wall surface and the plasma current). The dissipation time for these currents and the associated Joule dissipated energy should be properly accounted for to understand the energy balance of the plasma. It is not trivial because the resistivity of the conducting wall (liquid metal) is different from that of the plasma, the latter being subject to large changes as the temperature increases during the compression. This discussion is well beyond the scope of the present paper.

## Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## Acknowledgments

T. N. wishes to thank A. Tezcan for checking the analytical computations.

## Appendix A. Computation of $\oint d\theta R^{-2}$

$$ \begin{aligned} \oint \frac{\mathcal{J}}{R^2} d\theta &= \frac{\tilde{r}}{d\tilde{\psi}/d\tilde{r}} \frac{\varepsilon^2}{CB_0^\circ R_0^\circ} \oint \frac{1 + \varepsilon \tilde{\Delta}' \cos \theta}{1 + \varepsilon \tilde{r} \cos \theta + \varepsilon^2 \tilde{\Delta}} d\theta \\ &= \frac{\tilde{r}}{d\tilde{\psi}/d\tilde{r}} \frac{\varepsilon^2}{CB_0^\circ R_0^\circ} (1 + \varepsilon^2 \tilde{\Delta})^{-1} \oint \frac{1 + \alpha_1 \cos \theta}{1 + \alpha_2 \cos \theta} d\theta, \end{aligned} \quad (A.1) $$

with $\alpha_1 = \varepsilon \tilde{\Delta}'$ and $\alpha_2 = \varepsilon \tilde{r}/(1 + \varepsilon^2 \tilde{\Delta})$. The poloidal integral can be carried out using the residue theorem with the change of variable $z = e^{i\theta}$. This yields

$$ \oint \frac{1 + \alpha_1 \cos \theta}{1 + \alpha_2 \cos \theta} d\theta = - \frac{2i}{\alpha_2} \oint \frac{z + \alpha_1/2 (z^2 + 1)}{z(z - z_+) (z - z_-)} dz \quad (A.2) $$

with

$$ z_{\pm} = -\frac{1}{\alpha_2} \left( 1 \pm \sqrt{1 - \alpha_2^2} \right) \quad (A.3) $$

Note that $z_+z_- = 1$ and $z_- - z_+ = 2\sqrt{1 - \alpha_2^2}$. The poles that are within the unit circle are $z = z_-$ and $z = 0$. By the residue theorem, the integral is

$$ \oint \frac{1 + \alpha_1 \cos \theta}{1 + \alpha_2 \cos \theta} d\theta = 2\pi \left[ \frac{\alpha_1}{\alpha_2} + 2 \frac{z_- + \alpha_1/2 (z_-^2 + 1)}{\alpha_2 z_- (z_- - z_+)} \right] \quad (A.4) $$

The result can be plugged in any computer algebra software to yield a small aspect ratio expansion of the result. Using `sympy`, we find
$$ (1 + \varepsilon^2 \tilde{\Delta}) \left[ \oint \frac{1 + \alpha_1 \cos \theta}{1 + \alpha_2 \cos \theta} \right]^{-1} = \frac{1}{2\pi} \left[ 1 - \varepsilon^2 \left( \frac{\tilde{r}^2}{2} - \tilde{\Delta} - \frac{1}{2} \tilde{r} \tilde{\Delta}' \right) + \mathcal{O}(\varepsilon^4) \right] \quad (A.5) $$

## Appendix B. Neglect of $\psi_1$ term in the Shafranov shift equation

The reader who wants to check that the $\tilde{\psi}_1$ terms of (26) indeed lead to only $\varepsilon^2$ corrections to the flux solution of (24) can examine the right hand side of equation (27), had we kept them. For instance, the $-2\tilde{r}\tilde{\psi}'_0\tilde{\psi}'_1$ would yield, after replacing $\tilde{\psi}'_0$ with $\tilde{T}_0$, to the term $-2[\tilde{r}^2\tilde{T}_0/(f_0q_0)]\tilde{\psi}'_1$ on the right hand side of (27). The additional factor of $\varepsilon^2$ in the right hand side of equation (24) means that this would contribute to the coefficient of the $\tilde{\psi}'_1$ term to order $\varepsilon^2$, which we neglect. The term $-2[r\tilde{\psi}'_0\tilde{\psi}'_1\tilde{\Delta}'_0]'$ is of the same order, so that the use of equation (27) is justified.

## Appendix C. Comparison with Greene et al

The correction to the Furth-Yoshikawa scaling derived by Greene, Johnson and Weimer in reference [5] could, in principle, have been used to comment the results of refs. [11, 12, 13, 14], but it was not the case. Let us examine the equation proposed by ref. [5]. The pressure is expanded to second order in $\varepsilon$ and is written as $p = p^{(2)} + \varepsilon^2 p^{(4)}$. They use two distinct variables $\sigma \ge 1$ and $\tau \ge 1$ for the reduction of major radius and increase of the flux function $T$, respectively. To compare with our case, it is sufficient to use the prescription $\sigma = \tau = \sqrt{C}$. The comparison is made delicate because they have a different convention for the Shafranov shift. Instead of being a positive quantity which vanishes at the edge, it is a negative quantity which vanishes at the center. Also, the shift appearing in their expressions has been explicitly ordered with $\varepsilon^2$, so below we replace their symbol $\Delta$ with $\Delta/\varepsilon^2$. Then, their equation (113) for the pressure during the flux conserving compression reads.

$$ p^{(2)} = C^5 p^{(2)} \qquad (C.1) $$

$$ p^{(4)} = C^5 \left\{ p^{(4)} + \frac{5}{3} p^{(2)} \left[ 2 \frac{(C\Delta - \Delta^{\circ})}{\varepsilon^2 R_0^{\circ}} + r^{\circ} \frac{d\Delta/dr - d\Delta^{\circ}/dr}{\varepsilon^2 R_0^{\circ}} + g^{(2)} - g^{(2)} \right] \right\}, \quad (C.2) $$

where as usual the symbol $^{\circ}$ indicates quantities at the beginning of compression, $C = 1$, and the quantity we call $\tilde{T}$ is equal to $1 + \varepsilon^2 g^{(2)}$ in their notation. Therefore, above, we can replace $g^{(2)} - g^{(2)}$ with $C\tilde{T}_1/\varepsilon^2$ (confusing $C = 1$ with $C = 0$ as usual). Other substitutions are as follows:

$$ \frac{C\Delta - \Delta^{\circ}}{\varepsilon^2 R_0^{\circ}} = \tilde{\Delta}_1 \qquad (C.3) $$

$$ r^{\circ} \frac{d\Delta/dr - d\Delta^{\circ}/dr}{\varepsilon^2 R_0^{\circ}} = \tilde{r} \tilde{\Delta}'_1. \qquad (C.4) $$

The deviation with respect to the FY scaling comes from the square bracketed term in equation (C.2). We recognize the combination of Shafranov shift leading to the $f_1$ term in
our calculations. Let us denote it by $f_1^{\text{GJW}}$ to remember the difference in the convention taken for the Shafranov shift. This leads to the following expression, to be compared with (31):

$$ \left. \frac{\delta p}{Cp} \right|_{\text{GJW}} = \frac{5}{3} (2f_1^{\text{GJW}} + \tilde{T}_1). \qquad (C.5) $$

Actually, this expression can hardly be compared to ours. The problem is not with the sign difference so much as with the convention that $\Delta(r = 0) = 0$. This leads to no effect of the shift whatsoever on the magnetic axis, and one hardly gets physical insight from the $\tilde{T}_1$ term. This problem is due to the fact that if one writes the major radius, as Green et al. are doing:

$$ R^{\text{GJW}} = R_0^{\text{GJW}} - \varepsilon r \cos \theta - \varepsilon^2 \Delta^{\text{GJW}} \qquad (C.6) $$

with $\Delta^{\text{GJW}}(r = 0) = 0$, then it means that $R_0$ no longer represents the radius of the outermost flux surface, but the major radius of the magnetic axis. We don't have control over the latter quantity, since the profiles readjust according to Grad-Shafranov equation. Therefore, scaling the radius as $R^{\text{GJW}} = R^{\circ \text{GJW}}/C$ (see equation (106) in reference [5]) does not really make sense from an operational point of view. This explains the difference between our results and theirs.

Nonetheless, their expression suggests an alternative expression to equation (31). If we use equations (21)-(22) to express $\psi'_1/\psi'_0$, we get (at zero order in $\varepsilon$):

$$ \frac{\tilde{p}_1}{\tilde{p}_0} = \frac{5}{3} (-2f_1 + \tilde{T}_1) + \left( \frac{\tilde{p}'_0}{\tilde{p}_0} - \frac{q'_0}{q_0} \right) \frac{\tilde{\psi}_1}{\tilde{\psi}'_0} \qquad (C.7) $$

Up to the last term, proportional to $\tilde{\psi}_1$, this expression looks much like equation (C.5). Since $\tilde{\psi}_1(r = 0) = 0$, if we can neglect $\tilde{T}_1$, the result is that the net effect of the shift is twice that of the intuitive geometric effect described in the first paragraph of section 2.3. Examination of the zero shear case in section 3.3 and of the numerical results in sections 4.2 and 4.3 both support this interpretation.

## Appendix D. Axis shift

We can integrate (33) as follows:

$$
\begin{aligned}
f_1(0) &= \varepsilon^2 \tilde{\Delta}_1(0) \\
&= -\int_0^1 \varepsilon^2 \tilde{\Delta}'_1 \\
&= -\int_0^1 d\tilde{r} \frac{2q^2}{\tilde{r}^3} \int_0^{\tilde{r}} s^2 \tilde{p}'_0(s) ds \\
&= q^2 \int_0^1 d\tilde{r} \frac{d}{d\tilde{r}} \left( \frac{1}{\tilde{r}^2} \right) \int_0^{\tilde{r}} s^2 \tilde{p}'_0(s) ds \\
&= q^2 \left[ \int_0^1 \tilde{r}^2 \tilde{p}'_0(\tilde{r}) d\tilde{r} + \tilde{p}_0(0) \right] \\
&= q^2 \left[ \tilde{p}_0(0) - 2 \int_0^1 \tilde{r} \tilde{p}_0(\tilde{r}) d\tilde{r} \right] \\
&= \frac{q^2}{2} \left( \beta^\circ - \frac{\varepsilon^2}{q^2} \beta_P^\circ \right),
\end{aligned}
\qquad (D.1)
$$
where several integration by parts are carried out, and $\tilde{p}'_0(0) = 0$ is used to cancel one of the terms:

$$ \lim_{\tilde{r} \to 0} \frac{1}{\tilde{r}^2} \int_0^{\tilde{r}} s^2 \tilde{p}'_0(s) ds = 0. \qquad (D.2) $$

[1] Furth H P and Yoshikawa S 1970 *The Physics of Fluids* **13** 2593–2596 ISSN 0031-9171

[2] Laberge M 2019 *Journal of Fusion Energy* **38** 199–203 ISSN 0164-0313, 1572-9591

[3] Book D L, Cooper A L, Ford R, Gerber K A, Hammer D A, Jenkins D J, Robson A E and Turchi P J 1977 **3** 507–516

[4] Quimby D, Hoffman A and Vlases G 1981 *Nuclear Fusion* **21** 553

[5] Greene J M, Johnson J L and Weimer K E 1971 *The Physics of Fluids* **14** 671–683 ISSN 0031-9171

[6] Citrolo J C and Mullaney D H 1971 *IEEE Transactions on Nuclear Science* **18** 150–154 ISSN 1558-1578

[7] Bol K, Ellis R A, Eubank H, Furth H P, Jacobsen R A, Johnson L C, Mazzucato E, Stodiek W and Tolnas E L 1972 *Physical Review Letters* **29** 1495–1498

[8] Citrolo J C 1973 *IEEE Transactions on Nuclear Science* **20** 469–472 ISSN 1558-1578

[9] Shafranov V 1966 *Reviews of plasma physics (Edited by MA Leontovich)* vol. **2** p. 103

[10] Zakharov L and Shafranov V 1986 *Reviews of plasma physics (Edited by MA Leontovich)* vol. **11** p. 152

[11] Holmes J A, Peng Y K M and Lynch S J 1980 *Journal of Computational Physics* **36** 35–54 ISSN 0021-9991

[12] Holmes J A, Peng Y K M and Lynch S J 1980 *The Physics of Fluids* **23** 1874–1879 ISSN 0031-9171

[13] Albert D B 1980 *Nuclear Fusion* **20** 939–948 ISSN 0029-5515

[14] Brennan D, Froese A, Reynolds M, Barsky S, Wang Z, Delage M and Laberge M 2020 *Nuclear Fusion* **60** 046027

[15] Lütjens H, Bondeson A and Sauter O 1996 *Computer physics communications* **97** 219–260

[16] Chew G F, Goldberger M L, Low F E and Chandrasekhar S 1956 *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **236** 112–118

[17] Wesson J and Campbell D J 2011 *Tokamaks* vol 149 (Oxford university press)

[18] Blondel L, Sauter O and Merle A 2016 Axisymmetric tokamak equilibria computed with a predefined safety factor profile as chease input Tech. rep.

[19] Dory R A and Peng Y K M 1977 *Nuclear Fusion* **17** 21–31 ISSN 0029-5515

[20] Marinoni A, Sauter O and Coda S 2021 *Reviews of Modern Plasma Physics* **5** 6 ISSN 2367-3192

[21] Borowski S K and Kammash T 1979 *Energy* **4** 147–155 ISSN 0360-5442

[22] Clarke J F and Sigmar D J 1977 *Physical Review Letters* **38** 70–74

[23] Clarke J F and Sigmar D J 1976 Global properties of high pressure flux conserving tokamak equilibria Tech. Rep. ORNL/TM-5599 Oak Ridge National Lab., Tenn. (USA)

[24] Xue Ming-Lun 1981 A Combined Scheme of Adiabatic Compression and Accelerated Merging of Tokamak Tori