format: text

text:
# An improved unified solver for compressible and incompressible fluids involving free surfaces. II. Multi-time-step integration and applications

## 1 Introduction

This series of articles presents an improved solver for a challenging problem, the unified solution of compressible and incompressible fluids. After Harlow and Amsden proposed the ICE algorithm as a fully implicit solver for fluid equations in a conservative form , some approaches for this purpose have been investigated, such as the pressure-based semi-implicit algorithms , approaches based on the asymptotic expansion with respect to the local Mach number , and the boundary condition capturing method that treats compressible and incompressible materials separately . (See also recent reviews  for more details.) Among them, the CIP-CUP method, a pressure-based semi-implicit algorithm proposed by Yabe and Wang , has already been applied to many practical multi-material problems such as the laser machining of a metal plate , comet Shoemaker-Levy 9’s collision with the planet Jupiter , the interaction of a shock wave and a liquid drop , and the milk-crown formation on a liquid surface , and has been proven to be an efficient and robust solver for the unified solution . In this method, the Euler (or the Navier-Stokes) equations for fluid flows in a non-conservative form are selected as the governing equations, and the convection terms in the equations are solved explicitly by the CIP method , while the acoustic terms are solved implicitly by the CUP method . In Part I of this series , an improved solver for the convection terms was constructed using both an interpolation and an extrapolation function to describe the spatial profile of the density of materials. In the improved method, the cubic interpolation function used in the CIP scheme is replaced with the quadratic extrapolation functions constructed under some constraints for guaranteeing stability, only in a cell containing an interface between materials. This method allows us to solve the density interface with no dissipation across the interface, and is applicable to compressible flow problems, i.e., to variable-density problems.

As Part II of this series, this paper presents some further improvements for the CIP-CUP method, and also gives some application results. In Sec. 3, we attempt to incorporate the concept of the multi-time-step (MTS) integration techniques into the CIP-CUP algorithm. The MTS integration techniques, which allow us the efficient integration of dynamical systems involving some different time scales, have been used to solve a variety of problems, such as gravitational $`N`$-body problems , molecular dynamics , and atmospheric flow problems . This technique reduces the computational effort by integrating only rapidly varying components of a system (e.g., strongly accelerated particles or waves having fast phase speeds) with a small time interval, and others with a larger time interval. In the present study, based on this concept, different-nature terms in the fluid equations (convection, acoustic, and other terms) are solved using different time intervals. This modification allows us a more flexible determination of the time interval, and improves the accuracy and efficiency of the CIP-CUP. As will be shown, the adaptation of the MTS concept to the CIP-CUP integration is straightforward, because the CIP-CUP separately treats the terms having different natures, i.e., those having different characteristic time scales.

Section 4 presents a modified scheme to update the spatial derivatives during the computational steps except for the convection parts. As was shown in Part I, the spatial derivatives of the dependent variables are used explicitly as additional dependent variables to solve the convection terms. In the convection process, the derivatives are updated to obey the equations derived by taking derivatives of the governing equations . Those spatial derivatives should be updated by some way also in the non-convection processes. In Refs. , simple methods for this purpose have been proposed using classical centered differences. As was discussed in Part I, however, the use of such a classical discretization gives rise to the numerical dispersion and dissipation around the interfaces at which the spatial derivatives of the dependent variables tend to be discontinuous. We modify the centered scheme by adopting a simple extrapolation.

In Sec. 5, the applicability of the present methods to the compressible Navier-Stokes equations with a surface tension term is investigated, and in the Appendixes, some additional discussions regarding the averaging of the density (needed to solve the terms other than the convection) and the treatment of the surface-tension term are given. Section 6 presents some application results for the compressible and incompressible Kelvin-Helmholtz instabilities and the multibubble dynamics in an acoustic field, and Sec. 7 presents concluding remarks.

## 2 The CUP method and its variants

### 2.1 The CUP method

The CIP-CUP method  is a semi-implicit solver for the Euler equations of fluid flows:

$$\frac{\rho }{t}+𝐮\rho =\rho 𝐮,$$
(1)
$$\frac{𝐮}{t}+𝐮𝐮=\frac{p}{\rho }+\frac{𝐅}{\rho },$$
(2)
$$\frac{p}{t}+𝐮p=\rho C_S^2𝐮,$$
(3)
where $`\rho `$, $`𝐮`$, $`p`$, and $`C_S`$ denote the density, the velocity vector, the pressure, and the local sound speed, respectively, and F may contain the viscosity, the surface tension, and external forces. This method separately solves the terms in Eqs. (1)–(3) of different natures by a time splitting technique. The convection parts of these equations,

$$\frac{\rho }{t}+𝐮\rho =0,$$

$$\frac{𝐮}{t}+𝐮𝐮=0,$$

$$\frac{p}{t}+𝐮p=0,$$

are solved by the CIP method  (In the present work, of course, the hybrid interpolation-extrapolation method proposed in Part I is adapted), while the acoustic parts,

$$\frac{\rho }{t}=\rho 𝐮,$$
(4)
$$\frac{𝐮}{t}=\frac{p}{\rho },$$
(5)
$$\frac{p}{t}=\rho C_S^2𝐮,$$
(6)
are solved by the CUP method, which is an implicit finite difference method. The remaining part,

$$\frac{𝐮}{t}=\frac{𝐅}{\rho },$$
(7)
may be solved by some existing methods, such as the finite difference or finite volume method. (We refer to Eq. (7) as the “additional part”.) Solving these parts successively completes one step of the CIP-CUP time integration. In the following, the CUP method for the acoustic parts is concretely reviewed.

Discretizing the time derivatives on the LHS of Eqs. (4)–(6) and estimating the spatial derivatives on the RHS with future values, one obtains

$$\frac{\rho ^{n+1}\rho ^{}}{\mathrm{\Delta }t}=\rho ^{}𝐮^{n+1},$$
(8)
$$\frac{𝐮^{n+1}𝐮^{}}{\mathrm{\Delta }t}=\frac{p^{n+1}}{\rho ^{}},$$
(9)
$$\frac{p^{n+1}p^{}}{\mathrm{\Delta }t}=\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}𝐮^{n+1},$$
(10)
where $`n`$ is the number of the time step, $`\mathrm{\Delta }t`$ is the time interval, and the quantities with the superscript $``$ indicate the values after solving the parts other than the acoustic parts. (The acoustic parts should be solved at the final stage of a time step, because, as is shown below, solving these parts enforces the divergence-free condition for an incompressible fluid.) Taking divergence of Eq. (9) and substituting it into Eq. (10) yield the following pressure equation:

$$\frac{p^{n+1}p^{}}{\mathrm{\Delta }t}=\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}\mathrm{\Delta }t\frac{p^{n+1}}{\rho ^{}}\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}𝐮^{}.$$
(11)
(Substituting Eq. (6) into this and rewriting as

$$\left(\frac{p^{n+1}p^{}}{\mathrm{\Delta }t}\frac{p^{}}{t}\right)\mathrm{/}\mathrm{\Delta }t=\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}\frac{p^{n+1}}{\rho ^{}}$$

finds that Eq. (11) corresponds to a first-order approximation of the wave equation in terms of the pressure.) Almost the same pressure equation is given in Ref. , but discretization is done by a finite element technique.

After solving Eq. (11) to get $`p^{n+1}`$, the velocity is updated explicitly by Eq. (9). Also, the density is updated explicitly by

$$\rho ^{n+1}=\rho ^{}+\frac{p^{n+1}p^{}}{C_S^{\mathrm{\hspace{0.17em}2}}},$$
(12)
given by Eqs. (8) and (10). Generally, the spatial discretization for the above equations is performed using the 2nd-order centered finite differencing on the staggered grids.

When the sound speed is infinite, the pressure equation (11) is reduced to an elliptic equation,

$$\frac{p^{n+1}}{\rho ^{}}=\frac{𝐮^{}}{\mathrm{\Delta }t},$$
(13)
which corresponds to that used in the SMAC algorithm  for incompressible flows. This result reveals that the CUP method is applicable to both compressible and incompressible fluids.

In the case where the CIP or its variant is used to solve the convection parts, we also need to update the spatial derivatives, used explicitly in those schemes, in the non-convection processes. Simple schemes for this purpose have already been proposed . For a two-dimensional case, the scheme is represented as

$$_xf_{i,j}^{n+1}_xf_{i,j}^{}=\frac{f_{i+1,j}^{n+1}f_{i1,j}^{n+1}}{2h}\frac{f_{i+1,j}^{}f_{i1,j}^{}}{2h},$$
(14)
$$_yf_{i,j}^{n+1}_yf_{i,j}^{}=\frac{f_{i,j+1}^{n+1}f_{i,j1}^{n+1}}{2h}\frac{f_{i,j+1}^{}f_{i,j1}^{}}{2h},$$
(15)
where $`f`$ indicates an dependent variable, $`_xf=f/x`$, $`_yf=f/y`$, the subscripts $`i`$ and $`j`$ indicate $`(x,y)=(ih,jh)`$, and $`h`$ is the grid spacing (assumed to be constant for simplicity). In the case where the cross derivative $`_{xy}f`$ is used in the convection process , the following equation is additionally used:

$`_{xy}f_{i,j}^{n+1}_{xy}f_{i,j}^{}={\displaystyle \frac{f_{i+1,j+1}^{n+1}f_{i1,j+1}^{n+1}f_{i+1,j1}^{n+1}+f_{i1,j1}^{n+1}}{4h^2}}`$
$`{\displaystyle \frac{f_{i+1,j+1}^{}f_{i1,j+1}^{}f_{i+1,j1}^{}+f_{i1,j1}^{}}{4h^2}}.`$ (16)
Equations (14)–(2.1) can be solved explicitly using the quantities obtained before and after solving the non-convection parts.

### 2.2 Improved variants of the pressure equation

In 1994 , Ito proposed an improved variant of the pressure equation \[Eq. (11)\] by incorporating the concept of the exponential method for a heat-conduction equation . The 1-D formula of the variant is represented by

$$\frac{p^{n+1}p^{}}{\mathrm{\Delta }t}=\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}\mathrm{\Delta }t[\alpha \left(\frac{p_x^{n+1}}{\rho ^{}}\right)_x+(1\alpha )\left(\frac{p_x^{}}{\rho ^{}}\right)_x]\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}u_x^{},$$
(17)
where $`\alpha `$ is the weighting factor determined theoretically as

$$\alpha (E)=\frac{1}{1\mathrm{exp}(E)}\frac{1}{E},$$

$$E=\left(\frac{1}{\rho _{i+1/2}^{}}+\frac{1}{\rho _{i1/2}^{}}\right)\rho _i^{}\left(\frac{C_{S}^{}{}_{i}{}^{}\mathrm{\Delta }t}{h}\right)^2,$$

and $`\rho _{i\pm 1/2}^{}`$ is the density at $`x=h(i+1/2)`$, determined approximately with $`\rho _i^{}`$ and $`\rho _{i\pm 1}^{}`$ (see Appendix A). For $`E\mathrm{}`$ (i.e., $`C_{S}^{}{}_{i}{}^{}\mathrm{\Delta }t/h\mathrm{}`$), $`\alpha `$ converges to 1 and Eq. (17) is reduced to the elliptic equation (13); therefore, this variant is applicable to incompressible flows. For $`E0`$, $`\alpha `$ becomes $`1/2`$, resulting in a higher resolution than that of the conventional CUP . We adopted this variant in the present study.

The 2-D formula of this variant (not shown in Ito’s paper) may be represented by

$$\frac{p^{n+1}p^{}}{\mathrm{\Delta }t}=\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}\mathrm{\Delta }t(\stackrel{~}{G}1+\stackrel{~}{G}2)\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}𝐮^{},$$

where

$$\stackrel{~}{G}1=\alpha (E1)\left(\frac{p_x^{n+1}}{\rho ^{}}\right)_x+(1\alpha (E1))\left(\frac{p_x^{}}{\rho ^{}}\right)_x,$$

$$\stackrel{~}{G}2=\alpha (E2)\left(\frac{p_y^{n+1}}{\rho ^{}}\right)_y+(1\alpha (E2))\left(\frac{p_y^{}}{\rho ^{}}\right)_y,$$

$$E1=\left(\frac{1}{\rho _{i+1/2,j}^{}}+\frac{1}{\rho _{i1/2,j}^{}}\right)\rho _{i,j}^{}\left(\frac{C_{S}^{}{}_{i,j}{}^{}\mathrm{\Delta }t}{h}\right)^2,$$

$$E2=\left(\frac{1}{\rho _{i,j+1/2}^{}}+\frac{1}{\rho _{i,j1/2}^{}}\right)\rho _{i,j}^{}\left(\frac{C_{S}^{}{}_{i,j}{}^{}\mathrm{\Delta }t}{h}\right)^2.$$

We simplify this as follows to reduce the computational efforts:

$`{\displaystyle \frac{p^{n+1}p^{}}{\mathrm{\Delta }t}}=`$ $`\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}\mathrm{\Delta }t[\alpha (E_{\mathrm{max}}){\displaystyle \frac{p^{n+1}}{\rho ^{}}}+(1\alpha (E_{\mathrm{max}})){\displaystyle \frac{p^{}}{\rho ^{}}}]`$ (18)
$`\rho ^{}C_S^{\mathrm{\hspace{0.17em}2}}𝐮^{},`$
where

$$E_{\mathrm{max}}=\mathrm{max}(E1,E2).$$

This formula requires only one weighting factor. This simplification may be valid because $`E1E2`$ generally.

In Ref. , an alternative variant was proposed. Though this variant can provide higher resolution of the sound wave than that obtained by Ito’s variant, it is likely that the use of this variant as a part of the solver for fluid flows sometimes causes great violation of mass conservation . We therefore did not adapt it.

## 3 Multi-time-step integration

In many fluid flow problems, the time scale of a significant phenomenon is determined by the flow velocity (i.e., the characteristic speed of the convection parts), rather than by the sound speed or other characteristic speed. In certain cases, however, the acoustic parts (or others) have much smaller time scales (i.e., much greater characteristic speeds), resulting in an excessive restriction on the time interval, $`\mathrm{\Delta }t`$. The CIP-CUP overcomes this problem by solving the acoustic parts implicitly. Such an approach, however, guarantees only the stability, and probably provides inaccurate results in certain situations. For example, under an intermediate condition of flows such as weakly compressible flows, the acoustic parts may need to be solved by using a time interval of sufficiently small $`\kappa _C`$ that makes $`\kappa _u`$ very small, since the compressibility described by the acoustic parts plays an important role, where

$$\kappa _C\mathrm{max}(C_S)\mathrm{\Delta }t/h\mathrm{and}\kappa _u\mathrm{max}(\left|𝐮\right|)\mathrm{\Delta }t/h.$$

To achieve an efficient computation even in such a situation, we adopt the multi-time-step (MTS) integration technique to the CIP-CUP method.

The adaptation of the MTS integration to the CIP-CUP is straightforward because the components of different time scales (convection, sound propagation, and others) are treated separately. The time propagation performed in the CIP-CUP can be represented schematically as follows:

$$𝐔^{}=L1(\mathrm{\Delta }t)𝐔^n,$$
(19)
$$𝐔^{}=L2(\mathrm{\Delta }t)𝐔^{},$$
(20)
$$𝐔^{n+1}=L3(\mathrm{\Delta }t)𝐔^{},$$
(21)
where $`𝐔=(\rho ,𝐮,p)`$ and the operators $`L1(\mathrm{\Delta }t)`$, $`L2(\mathrm{\Delta }t)`$, and $`L3(\mathrm{\Delta }t)`$ indicate the discrete propagators for the corresponding integration steps (the convection, the additional, and the acoustic steps, respectively). These equations can be summarized as

$$𝐔^{n+1}=L3(\mathrm{\Delta }t)L2(\mathrm{\Delta }t)L1(\mathrm{\Delta }t)𝐔^n.$$
(22)
If we factor $`L2(\mathrm{\Delta }t)`$ and $`L3(\mathrm{\Delta }t)`$ into some identical pieces, we get

$$𝐔^{n+1}=[L3(\mathrm{\Delta }t/m3)]^{m3}[L2(\mathrm{\Delta }t/m2)]^{m2}L1(\mathrm{\Delta }t)𝐔^n,$$
(23)
where $`m2`$ and $`m3`$ are positive integers. This equation means that, after solving the convection parts with $`\mathrm{\Delta }t`$, the additional part is solved $`m2`$ times with $`\mathrm{\Delta }t/m2`$, and subsequently the acoustic parts are solved $`m3`$ times with $`\mathrm{\Delta }t/m3`$. If necessary, we divide $`L2`$ into two parts as $`L2=L2^{(2)}L2^{(1)}`$, and modify Eq. (23) as

$$𝐔^{n+1}=[L3(\mathrm{\Delta }t/m3)L2^{(2)}(\mathrm{\Delta }t/m3)]^{m3}[L2^{(1)}(\mathrm{\Delta }t/m2)]^{m2}L1(\mathrm{\Delta }t)𝐔^n$$
(24)
or

$`𝐔^{n+1}=[L3(\mathrm{\Delta }t/m3)L2^{(2)}(\mathrm{\Delta }t/m3)[L2^{(1)}(\mathrm{\Delta }t/(m2m3))]^{m2}]^{m3}L1(\mathrm{\Delta }t)𝐔^n.`$
Such a treatment of the additional parts is necessary, e.g., when the surface-tension term exists. (The surface-tension part, represented by $`L2^{(2)}`$, should be solved together with the acoustic parts because the surface tension always needs to balance with the pressure jump at the interface.) Some other forms of factorization adjusted for a governing system can be employed.

We determine the fundamental time interval, $`\mathrm{\Delta }t`$, by

$$\mathrm{\Delta }t=\mathrm{min}(c1\frac{h}{\mathrm{max}(\left|𝐮\right|)},c2\frac{h}{\mathrm{max}(C_S)}),$$

where the parameters are typically set to $`c1=0.2`$ and $`c2=10`$. $`m3`$ is determined so that, e.g., the following condition is satisfied:

$$\frac{\mathrm{\Delta }t}{m3}\mathrm{min}(c3\frac{h}{\mathrm{max}(\left|C_{S}^{}{}_{i+1,j}{}^{}C_{S}^{}{}_{i,j}{}^{}\right|)},c3\frac{h}{\mathrm{max}(\left|C_{S}^{}{}_{i,j+1}{}^{}C_{S}^{}{}_{i,j}{}^{}\right|)},\mathrm{\Delta }t),$$

where the spatial difference of sound speed is used as a criterion. If the surface-tension term exists, the stability condition needed to solve it ,

$$\frac{\mathrm{\Delta }t}{m3}\mathrm{\Delta }t_{st}\sqrt{\frac{(\rho _a+\rho _b)h^3}{4\pi \sigma }},$$
(26)
should be taken into consideration to determine $`m3`$, where $`\rho _a`$ and $`\rho _b`$ are the densities of materials on different sides of the interface, and $`\sigma `$ is the surface-tension coefficient. Also, $`m2`$ is determined based on the stability condition of a method for the viscous term or others.

## 4 Stabilization of the derivative advancement

At the interfaces where the value or derivatives of physical quantities are in general discontinuous, the use of the conventional centered schemes (14)–(2.1) for updating the spatial derivatives is not suitable, providing less accurate results, as will be demonstrated below using a numerical example. To modify the schemes, we recall here the extrapolation concept investigated in Part I  and by others . As has been pointed out in Ref. , an interpolation or differencing across the phase boundary may not physically valid, and gives rise to serious numerical errors. In the following, those schemes are modified by adapting a simple extrapolation.

Equation (14) can be rewritten as

$$_xf_{i,j}^{n+1}_xf_{i,j}^{}=\frac{1}{2}\left(\frac{d_{i+1,j}d_{i,j}}{h}+\frac{d_{i,j}d_{i1,j}}{h}\right),$$
(27)
$$d_{i,j}f_{i,j}^{n+1}f_{i,j}^{}.$$

We introduce here a switching parameter $`H`$, defined as

$$H_{i+1/2,j}=\{\begin{array}{cccccccccccccccccccc}1& \mathrm{for}\varphi _{i+1,j}^{n+1}\varphi _{i,j}^{n+1}>0,& & & & & & & & & & & & & & & & & & \\ 0& \mathrm{otherwise},& & & & & & & & & & & & & & & & & & \end{array}$$
(28)
$$H_{i,j+1/2}=\{\begin{array}{cccccccccccccccccccc}1& \mathrm{for}\varphi _{i,j+1}^{n+1}\varphi _{i,j}^{n+1}>0,& & & & & & & & & & & & & & & & & & \\ 0& \mathrm{otherwise},& & & & & & & & & & & & & & & & & & \end{array}$$
(29)
where $`\varphi `$ is the ID function (a density function or a level set function) defined as $`\varphi >0`$ in the region occupied by a material and $`\varphi <0`$ elsewhere, and is updated by solving

$$\frac{\varphi }{t}+𝐮\varphi =0.$$
(30)
Using this, we modify Eq. (27) as

$$_xf_{i,j}^{n+1}_xf_{i,j}^{}=\frac{1}{2}\left(H_{i+1/2,j}\frac{d_{i+1,j}d_{i,j}}{h}+H_{i1/2,j}\frac{d_{i,j}d_{i1,j}}{h}\right).$$
(31)
This modification vanishes the derivative $`(d_{i+1,j}d_{i,j})/h`$ when $`\varphi _{i+1,j}^{n+1}\varphi _{i,j}^{n+1}<0`$, resulting in the advancement of $`_xf_{i,j}`$ using only the values of an identical material. This procedure corresponds to an extrapolation such as $`d_{i+1,j}=d_{i,j}`$. Equation (15) can be modified in the same manner, as

$$_yf_{i,j}^{n+1}_yf_{i,j}^{}=\frac{1}{2}\left(H_{i,j+1/2}\frac{d_{i,j+1}d_{i,j}}{h}+H_{i,j1/2}\frac{d_{i,j}d_{i,j1}}{h}\right).$$
(32)
The scheme (2.1) for the cross derivative can be rewritten as

$`_{xy}f_{i,j}^{n+1}_{xy}f_{i,j}^{}={\displaystyle \frac{1}{4h^2}}[(d_{i+1,j+1}d_{i,j})(d_{i1,j+1}d_{i,j})`$
$`(d_{i+1,j1}d_{i,j})+(d_{i1,j1}d_{i,j})].`$ (33)
We modify this as

$`_{xy}f_{i,j}^{n+1}`$ $`_{xy}f_{i,j}^{}={\displaystyle \frac{1}{4h^2}}`$ (34)
$`\times [H_{i+1/2,j+1/2}(d_{i+1,j+1}d_{i,j})H_{i1/2,j+1/2}(d_{i1,j+1}d_{i,j})`$
$`H_{i+1/2,j1/2}(d_{i+1,j1}d_{i,j})+H_{i1/2,j1/2}(d_{i1,j1}d_{i,j})],`$
where

$$H_{i\pm 1/2,j\pm 1/2}=\{\begin{array}{cccccccccccccccccccc}1& \mathrm{for}\varphi _{i\pm 1,j\pm 1}^{n+1}\varphi _{i,j}^{n+1}>0,& & & & & & & & & & & & & & & & & & \\ 0& \mathrm{otherwise}.& & & & & & & & & & & & & & & & & & \end{array}$$
(35)
The extrapolations used above are merely rough compared with those used to solve the convection parts ; these, however, might be sufficient, because the derivatives require lower accuracy than those required for the non-derivatives represented here by $`f`$.

Let us perform a numerical test to demonstrate the effectiveness of the above modification. A one-dimensional nonlinear sound propagation is solved based on Eqs. (4)–(6). The convection and other components are neglected for simplicity. The initial condition is

$$\begin{array}{c}\rho (x,0)=\{\begin{array}{cccccccccccccccccccc}1.025\times 10^3& \mathrm{for}x<0.5,& & & & & & & & & & & & & & & & & & \\ 1& \mathrm{elsewhere},& & & & & & & & & & & & & & & & & & \end{array}\hfill \\ u(x,0)=0,p(x,0)=1,\hfill \\ x[0,\mathrm{}),\hfill \end{array}$$

and the boundary condition is

$$p(0,t)=1.10.1\mathrm{cos}\omega _1t\mathrm{and}u(0,t)/x=0,$$

where the angular frequency $`\omega _1`$ is determined so that the wavelength of the emitted sound wave is $`0.225`$. The initial values of all derivatives are zero. The square of sound speed is determined by

$$C_S^2=\{\begin{array}{cccccccccccccccccccc}7(p+3172.04)/\rho & \mathrm{for}x<0.5,& & & & & & & & & & & & & & & & & & \\ 1.4p/\rho & \mathrm{elsewhere}.& & & & & & & & & & & & & & & & & & \end{array}$$

Other parameters are $`\mathrm{\Delta }t=2\times 10^5`$ and $`h=5\times 10^3`$. (Under this condition, $`\kappa _C`$ is about 0.6 for $`x0.5`$ and is about 0.15 elsewhere.) The ID function is set to be a color function,

$$\varphi =\{\begin{array}{cccccccccccccccccccc}1& \mathrm{for}x<0.5,& & & & & & & & & & & & & & & & & & \\ 1& \mathrm{elsewhere}.& & & & & & & & & & & & & & & & & & \end{array}$$

The matrix equation given by discretizing Eq. (18) is solved by the Red-Black Gauss-Seidel method with the convergence criterion of $`\epsilon _{res}<10^6`$, where

$$\epsilon _{res}=p_i^{n+1,(k)}p_i^{n+1,(k1)}_2\mathrm{/}p_i^{}_2,$$

and $`k`$ denotes the number of the Gauss-Seidel iteration. Figures 1(a) and 1(b) show the calculated pressure and density distributions at $`t=1.602\times 10^2`$. Also, Figs. 1(c) and 1(d) show the density gradients at the same $`t`$, calculated using the conventional and the modified schemes. (In the acoustic step, the derivatives are not used to update the non-derivatives; thus, the results shown in Figs. 1(a) and 1(b) are irrespective of the schemes for the derivatives.) As can be clearly seen, the result obtained by the conventional scheme has a strong overshoot at the interface, while the result obtained by the modified scheme is smooth. The result given using the modified scheme with $`\mathrm{\Delta }t/4`$ and $`h/4`$ (the dashed lines) is similar to the latter result. These results can roughly prove the validity of the modified scheme, although more detailed discussions should be necessary.

## 5 Application to the compressible Navier-Stokes equations

In this section, as a summary of methods, we illustrate how to apply the methods discussed in this series to the compressible Navier-Stokes equations with a surface tension term:

$`{\displaystyle \frac{\rho }{t}}+𝐮\rho `$ $`=`$ $`\rho 𝐮,`$ (36)
$`{\displaystyle \frac{𝐮}{t}}+𝐮𝐮`$ $`=`$ $`{\displaystyle \frac{p}{\rho }}+{\displaystyle \frac{1}{\rho }}\left(2(\mu 𝐓){\displaystyle \frac{2}{3}}(\mu 𝐮)\right)+{\displaystyle \frac{𝐅_{st}}{\rho }},`$ (37)
$`{\displaystyle \frac{p}{t}}+𝐮p`$ $`=`$ $`\rho C_s^2𝐮,`$ (38)
where $`𝐓`$ is the deformation tensor, $`\mu `$ is the viscosity coefficient, and $`𝐅_{st}`$ indicates the surface tension as a volume force. We factor them into the following four systems:

Convection system:

$$\{\frac{\rho }{t}+𝐮\rho =0,\frac{𝐮}{t}+𝐮𝐮=0,\frac{p}{t}+𝐮p=0,$$

Viscous system:

$$\{\frac{\rho }{t}=0,\frac{𝐮}{t}=\frac{1}{\rho }(2(\mu 𝐓)\frac{2}{3}(\mu 𝐮)),\frac{p}{t}=0,$$

Surface-tension system:

$$\{\frac{\rho }{t}=0,\frac{𝐮}{t}=\frac{𝐅_{st}}{\rho },\frac{p}{t}=0,$$

Acoustic system:

$$\{\frac{\rho }{t}=\rho 𝐮,\frac{𝐮}{t}=\frac{p}{\rho },\frac{p}{t}=\rho C_s^2𝐮.$$

Successively solving these four systems by the MTS technique completes one step of the time integration. Among the convection equations, one in terms of the density is solved by the hybrid interpolation-extrapolation method introduced in Part I, and the others are solved by the conventional CIP (the “Type B” multidimensional formula  is adapted). The time integration of $`\varphi `$ is performed together with the convection system. The viscous system is solved by the conventional second-order centered finite differencing, and the surface-tension system is solved by the CSF model , coupled with the smoothing procedure introduced in Appendix B. The acoustic system is solved by the modified CUP method reviewed in Sec. 2.2.

## 6 Application results

This section presents the application results for the multiphase flow problems, given using the present methods. Compressible and incompressible Kelvin-Helmholtz instabilities and multibubble dynamics in an acoustic field were selected as application examples.

Compressible and incompressible Kelvin-Helmholtz instabilities. It is well known that an interface between fluids of different densities is unstable. When the velocity field around the interface is perturbed, the interface forms a complicated structure . This instability, called the Kelvin-Helmholtz instability or the Rayleigh-Taylor instability, has been the subject of quite a number of theoretical, experimental, and numerical studies , and even the CIP have been employed to investigate it . Using this example, the effectiveness of the extrapolation technique for the convection parts and that of the MTS technique are demonstrated. The initial arrangement of materials is shown in Fig. 2. The density of the heavy material (one for $`x<0`$) is $`\rho _h=1`$, and that of the light one $`\rho _l=1/3`$, and the initial pressure is $`p_0(x,y,t=0)=1`$. The initial velocity field \[$`𝐮(x,y,t=0)=(u_0(x,y),v_0(x,y))`$\] is set to

$$u_0(x,y)=sgn(x)u_{\mathrm{max}}\mathrm{exp}(2\pi \left|x\right|/\lambda )\mathrm{sin}(2\pi y/\lambda ),$$

$$v_0(x,y)=u_{\mathrm{max}}\mathrm{exp}(2\pi \left|x\right|/\lambda )\mathrm{cos}(2\pi y/\lambda ),$$

$$(x,y)[1.25,1.25]\times [0,0.5],$$

which satisfies $`𝐮=0`$, where $`u_{\mathrm{max}}`$ is the maximum velocity and $`\lambda `$ ($`=1`$) is the wavelength of the perturbation. The non-dimensional parameter $`\lambda \sqrt{\rho _hp_0}/\mu `$ is fixed to 1400, and $`\mu `$ is assumed to be constant. The surface tension and the gravity are neglected. The slip boundary condition is adapted to the boundaries of $`y=0`$ and $`y=0.5`$, and the Neumann boundary condition ($`/\stackrel{}{n}=0`$) to $`x=1.25`$ and $`x=1.25`$. Figures 3(a)–3(d) show the time sequence of $`\rho `$ for $`u_{\mathrm{max}}=0.7`$, $`C_S^2=1.4p/\rho `$ (the Mach number under this condition is about 0.6), and $`h`$ $`(=\mathrm{\Delta }x=\mathrm{\Delta }y)=0.5/60`$, at $`t=4`$. Here, the linear monotone function  is used for the extrapolation, and Eq. (23) is used for the MTS integration. The parameters for the time intervals are $`c1=0.5`$, $`c2=\mathrm{}`$, and $`c3=0.1`$. (The inner iteration number for the acoustic parts, $`m3`$, is $`713`$ in this case.) These results present sharp descriptions of the density interface, while the numerical diffusion and oscillation can be observed in the result given when the convection parts are solved by the conventional CIP (Figs. 3(e) and 4). Figures 5(a)–5(d) show $`\varphi =0`$ surfaces at $`t=4`$ for different $`c3`$. (For $`c3=\mathrm{}`$, the MTS integration is disabled, that is, a single time step is used.) As is clearly seen, decreasing $`c3`$ changes the result. A similar result to that for $`c3=0.1`$ can be given, even using a smaller $`\mathrm{\Delta }t`$ (see Figs. 5(e) and 5(f), which show the results for $`c1=0.2`$ and different $`c3`$). These results reveal that integrating the acoustic parts with a smaller time interval improves the accuracy of the solution, even if the time interval for the convection parts holds.

While the above example can be solved by an explicit scheme without any difficulty, the next example is problematical. Figure 6 shows the results given by artificially resetting the sound speed to $`C_S^2=1000\times 1.4p/\rho `$ for only the heavy material. This setting results in a much greater sound speed than the flow velocity, and also results in the coexistence of materials of greatly different compressibilities. The parameters for the time integration are $`c1=0.25`$, $`c2=\mathrm{}`$, and $`c3=5.0`$. (Under this setting, $`m3=35`$ and $`\kappa _C/m345`$.) This result is not changed noticeably by decreasing the time interval, proving the applicability of the present method to a high $`\kappa _C`$ condition and proving its robustness. Using this example, we perform here a convergency test. Figures 7(a)–7(c) show the interfaces at $`t=4`$ given using the same parameters except for the grid width set to $`h=h_0`$, $`h_0/1.5`$, or $`h_0/2`$, respectively, where $`h_0=0.5/60`$. The refined results are in good agreement with that for $`h=h_0`$. In contrast, the result given by solving the convection parts with only the interpolation (i.e. by the conventional CIP) changes as the computational grids are refined (see Figs. 7(d)–7(f)); the convergency of the solution is obviously lower than that obtained by the hybrid method. These results show that the improvement given in Part I accelerates the convergency, and that the hybrid scheme can accurately describe such deformable interfaces.

Two-bubble dynamics in an acoustic field. When a sound wave is applied, a bubble immersed in a liquid begins volume oscillation. When other bubbles exist, they interact acoustically with each other, resulting in the changes in the oscillation amplitude and phase and the effective resonance frequencies . Moreover, it is known that an acoustic interaction force called the secondary Bjerknes force acts between such pulsating bubbles . The force is attractive when the bubbles pulsate in-phase, and is repulsive otherwise. Such effects resulting from the radiative interaction between bubbles have been important subjects not only for physicists and mechanical and acoustical engineers, but also for medical engineers and chemists , and novel insights regarding these effects have been provided even in very recent years . (In Refs. , for example, the author discovered that in multibubble cases, the phase shift of the bubbles’ pulsation can take place not only at their natural frequencies, but also at some other driving frequencies. In a more recent paper , it was found that the latter characteristic frequencies cause the sign reversal of the secondary Bjerknes force .)

Solving this problem numerically has some interesting points. (1) The gas inside the bubbles and the liquid surrounding them have quite different densities and compressibilities. (2) The compressibility of the gas is not negligible in principle, meaning that a solver for incompressible flows is not applicable. (3) Especially for small bubbles, the viscosity of the surrounding liquid and the surface tension are not negligible. (4) Many different time scales exist, determined by the flow velocity, the sound speed, the surface tension, the viscosity, the bubbles’ natural frequencies, and the frequency of an external sound. (5) The topology of the interfaces changes when the bubbles coalesce as a result of their radiative interaction. Based on these interesting factors, we consider this problem to be a very good example for demonstrating the abilities of our methods.

The axisymmetric coordinate $`(r,z)`$ is selected for the computational domain, and the centers of the bubbles are located on the central axis. Bubbles 1 and 2, whose equilibrium radii are $`R_{10}=5`$ $`\mu `$m and $`R_{20}=9`$ $`\mu `$m, respectively, are filled with a gas of the constant specific heat ratio $`\gamma =1.33`$. The equilibrium density and the viscosity coefficient are 1.226 kg/m<sup>3</sup> and $`1.78\times 10^5`$ Pa s for the gas phase, and 1000 kg/m<sup>3</sup> and $`1.137\times 10^3`$ Pa s for the liquid phase. The surface tension coefficient is $`\sigma =7.28\times 10^2`$ Pa m, and the equilibrium pressure of the liquid is $`P_0=101.3`$ kPa. The initial pressures inside the bubbles are thus determined by $`P_0+2\sigma /R_{j0}`$ (for $`j=1`$ and $`2`$), in order to balance with the surface tension, where $`2\sigma /R_{10}0.29P_0`$ and $`2\sigma /R_{20}0.16P_0`$. The sound speed is set to $`C_S^2=\gamma p/\rho `$ for the gas, and to $`C_S^2=7(p+3172.04P_0)/\rho `$ for the liquid. The external sound is applied as the boundary condition to the pressure by $`p(\mathrm{\Gamma })=P_0(1+0.3\mathrm{sin}\omega t)`$, where $`\mathrm{\Gamma }`$ denotes the boundaries of the computational domain except for $`r=0`$, and $`\omega `$ is the angular frequency of the sound. The boundary condition for the velocity is free. The initial distance between the centers of the bubbles \[$`D(t=0)`$\] is fixed to $`20`$ $`\mu `$m. The grid width assumed to be constant is set to $`h=(5/20)`$ $`\mu `$m, and the total number of grids is $`100\times 310`$. $`\rho _a+\rho _b`$ in Eq. (26) is set to 1000 kg/m<sup>3</sup> (= the equilibrium density of the liquid), because the density of the gas is negligible, and that of the liquid can be assumed to be almost constant in time.

Figure 8 shows the bubbles’ mean radii and mass centers calculated for $`\omega =\omega _{10}`$ and $`\omega =\omega _{10}/1.8`$ ($`\omega _{20}`$), where $`\omega _{j0}`$ is the linear, monopole natural frequency of bubble $`j`$ determined by

$$\omega _{j0}=\sqrt{\left[3\gamma P_0+(3\gamma 1)2\sigma /R_{j0}\right]/\rho R_{j0}^2},$$

and the parameters for the MTS integration were set to $`c1=0.2`$, $`c2=20`$, and $`c3=2`$, resulting in $`m38`$. Also, Figs. 9 and 10 display the interfaces at selected times. We can observe the repulsion of the bubbles when $`\omega =\omega _{10}`$, while the attraction when $`\omega =\omega _{10}/1.8`$. These results are validated by the theory presented in Ref. , and reexamined in Refs. . The theory determines the sign of the interaction force ($`𝐅_{2B}`$) by

$$sgn(𝐅_{2B})=sgn(\mathrm{cos}(\varphi _1\varphi _2)),$$
(39)
where

$$\varphi _1=\mathrm{tan}^1\left(B_1/A_1\right)[0,2\pi ]$$

with

$$A_1=\frac{H_1F+M_2G}{F^2+G^2},B_1=\frac{H_1GM_2F}{F^2+G^2},$$

$$F=L_1L_2\frac{R_{10}R_{20}}{D^2}\omega ^4M_1M_2,$$

$$G=L_1M_2+L_2M_1,H_1=L_2+\frac{R_{20}}{D}\omega ^2,$$

$$L_1=(\omega _{10}^2\omega ^2),L_2=(\omega _{20}^2\omega ^2),$$

$$M_1=\delta _1\omega ,M_2=\delta _2\omega ,$$

where $`\delta _j`$, set to $`(4\mu _l/\rho _lR_{j0}^2)+(\omega ^2R_{j0}/c_l)`$ in the present study, is the damping coefficient , $`\mu _l`$, $`\rho _l`$, and $`c_l`$ are the viscosity coefficient, the equilibrium density and the sound speed, respectively, of the liquid, and exchanging 1 and 2 (or 10 and 20) in the subscripts of these equations yields the expression for $`\varphi _2`$. The positive sign of Eq. (39) indicates attraction, while the negative indicates repulsion. Figure 11 shows $`sgn(𝐅_{2B})`$ as a function of $`\omega `$ for $`D=20`$ $`\mu `$m, revealing that the bubbles repel each other when $`\omega =\omega _{10}`$, while they attract when $`\omega =\omega _{10}/1.8`$. This theoretical prediction is well reproduced in the numerical results. (More detailed discussions from the physical viewpoint are given in Ref. .) Meanwhile, in the result for $`\omega =\omega _{10}/1.8`$, we can observe the coalescence of the bubbles as a result of the attraction. This result proves the ability of the present methods to treat the topology change.

Figure 12 shows $`\mathrm{\Delta }t`$ normalized by $`c1h/\mathrm{max}(\left|𝐮\right|)`$, $`c2h/\mathrm{max}(C_S)`$, or $`\mathrm{\Delta }t_{st}`$, as functions of time. This figure reveals that the fundamental time interval was mainly determined by the sound speed, whereas the flow velocity was responsible when a rapid motion of the interfaces, caused by the coalescence, occurred.

Using the example for $`\omega =\omega _{10}`$, we investigate here the effectiveness of the MTS integration on this problem. Figure 13 shows the bubbles’ mean radii as functions of time, for $`c1=0.2`$ and different $`c2`$ and $`c3`$. Large differences cannot be seen between the results for $`c34`$, while the results for $`c3=8`$ ($`m32`$) are obviously inaccurate. This result proves that the MTS integration has obviously contributed to the accurate solutions given for $`c34`$.

During the computation for $`\omega =\omega _{10}/1.8`$ introduced above, we sometimes observed negative $`p`$ near the interfaces. In such a case, we adopted $`p=\mathrm{max}(p,\epsilon _cP_0)`$, where $`\epsilon _c`$ is set to $`10^3`$; the numerical results were insensitive to the magnitude of this parameter. This problem would be overcome by employing the oscillation-free variants of the CIP .

## 7 Conclusions

In this series of articles, we have proposed an improved unified solver for compressible and incompressible fluids involving free surfaces, based on the CIP-CUP method, by adapting several improvements and modifications. (The most significant one given here in Part II is the adaptation of the MTS integration technique, which makes the determination of the time interval very flexible.) High accuracy and excellent robustness of the improved methods have been demonstrated by using examples of free-surface flows that contain both compressible and incompressible materials. The present methods, however, face the following challenges:

1. Optimizing the extrapolation function for the convection parts

In Part I of this series, we proposed five kinds of extrapolation functions used to solve the convection of the density. Although a concrete discussion has not been provided in the present paper, we observed that the most accurate result for a different problem was achieved by a different extrapolation function. The optimization of the function based on some criterion is, therefore, sorely expected.

2. Optimizing the time intervals

Although in the present study we have empirically determined the parameters for the MTS integration, an automatic, optimized determination should be useful in a practical application. The optimal time interval for the acoustic parts might be determined by a criterion based on the compressibility of fluids. If, for example, materials can be considered to be almost completely incompressible, as has been well known previously, one can use a time interval of an infinite CFL number with respect to the sound speed, whereas when a compressible material exists, as has been demonstrated in this paper, one needs a time interval depending on the sound speed. This means that the maximum values of $`c2`$ and $`c3`$ sufficient for accurate computations are rightly dependent on the compressibility.

The above subjects will be addressed in a future paper.

Acknowledgement

The author would like to thank Dr. N. Masuda at Gunma University for his helpful supply of materials on the multi-time-step integration methods (or the individual timestep schemes) used in the astrophysical community.

APPENDIX A: Averaging at phase boundary

To discretize the acoustic, viscous, and surface tension terms on the staggered grids, we need the values of $`\rho `$ at the velocity positions, i.e., $`\rho _{i+1/2,j}`$ and $`\rho _{i,j+1/2}`$. In this appendix, we briefly discuss how to estimate them.

As discussed in Part I, the phase boundary is recognized using the zero level set of the ID function (Identification Function; the color or the level set function), $`\varphi `$, and materials are identified using the sign of the function. If $`\varphi _{i+1,j}\varphi _{i,j}<0`$ is true, it is recognized that the call between ($`x_{i+1,j},y_{i+1,j}`$) and ($`x_{i,j},y_{i,j}`$) contains an interface. In such a cell, we estimate the densities by the following VOF -like procedure, which implies the weighted linear interpolation:

$$\rho _{i+1/2,j}^{}=\frac{\left|\varphi _{i+1,j}^{n+1}\right|\rho _{i+1,j}^{}+\left|\varphi _{i,j}^{n+1}\right|\rho _{i,j}^{}}{\left|\varphi _{i+1,j}^{n+1}\right|+\left|\varphi _{i,j}^{n+1}\right|}\mathrm{for}\varphi _{i+1,j}^{n+1}\varphi _{i,j}^{n+1}<0,$$
(40)
$$\rho _{i,j+1/2}^{}=\frac{\left|\varphi _{i,j+1}^{n+1}\right|\rho _{i,j+1}^{}+\left|\varphi _{i,j}^{n+1}\right|\rho _{i,j}^{}}{\left|\varphi _{i,j+1}^{n+1}\right|+\left|\varphi _{i,j}^{n+1}\right|}\mathrm{for}\varphi _{i,j+1}^{n+1}\varphi _{i,j}^{n+1}<0,$$
(41)
where $`\varphi ^{n+1}`$ ($`=\varphi ^{}=\varphi ^{}`$. Note that the position of the interface is changed only when solving the convection parts.) is the ID function after solving the convection part. In the remaining cells, we use the simple average,

$$\rho _{i+1/2,j}^{}=\frac{\rho _{i+1,j}^{}+\rho _{i,j}^{}}{2},$$

$$\rho _{i,j+1/2}^{}=\frac{\rho _{i,j+1}^{}+\rho _{i,j}^{}}{2},$$

which has been used in the conventional CIP algorithm , or in others.

The above averaging scheme is also used to estimate the values of some other quantities at the velocity positions, such as the viscosity coefficient.

APPENDIX B: On the surface tension term

This appendix concerns how to calculate the surface tension term in our code.

In the model named the CSF (continuum surface force), proposed by Brackbill et al , the surface tension as a volume force is represented by

$$𝐅_{st}=\sigma \kappa \theta ,$$
(42)
where $`\sigma `$ is the surface tension coefficient, $`\kappa `$ defined as

$$\kappa =\left(\frac{\theta }{\left|\theta \right|}\right)$$

is the mean radius of the surface curvature, and $`\theta `$ denotes the color function defined as $`\theta =1`$ for a material, and $`\theta =0`$ otherwise. The singularity appearing in $`\theta `$ at the interface is mollified by some smoothing procedures . In the present study, the smoothed $`\theta `$ is constructed as follows:

Step 1: The non-smoothed color function is made from the ID function, which is an arbitrary function whose $`\varphi =0`$ surface represents the interface , by setting $`\theta _{i,j}=1`$ for $`\varphi _{i,j}>0`$ and $`\theta _{i,j}=0`$ for $`\varphi _{i,j}<0`$. ($`\varphi `$ always has a non-zero value at the gird points .)

Step 2: The values of the color function at the velocity positions, $`\theta _{i+1/2,j}`$ and $`\theta _{i,j+1/2}`$, are determined by the averaging procedure described in Appendix A.

Step 3: Using the 4-point simple average, the weakly smoothed values of the color function at the original position ($`i,j`$) is obtained as

$$\overline{\theta }_{i,j}=0.25(\theta _{i+1/2,j}+\theta _{i1/2,j}+\theta _{i,j+1/2}+\theta _{i,j1/2}).$$

Step 4: Using the following 5-point smoother, $`\overline{\theta }`$ is smoothed further:

$$\overline{\theta }_{i,j}^{(m+1)}=(1\epsilon )\overline{\theta }_{i,j}^{(m)}+\epsilon \frac{1}{4}(\overline{\theta }_{i+1,j}^{(m)}+\overline{\theta }_{i1,j}^{(m)}+\overline{\theta }_{i,j+1}^{(m)}+\overline{\theta }_{i,j1}^{(m)}),$$

where $`m`$ ($`=0,1,\mathrm{},M1`$) is the number of iteration, $`M`$ is the total number of iteration, and $`\epsilon `$ is a small positive value.

Steps 2 and 3 are necessary to approximately take into account the phase property of the interface, which is not contained in $`\theta `$ given at Step 1, but is latent between the grid points. $`\epsilon `$ and $`M`$ are typically set to $`\epsilon =0.15`$ and $`M=20`$.

Replacing $`\theta `$ with $`\overline{\theta }^{(M)}`$, we calculate the surface-tension term by the second-order centered finite differencing on the staggered grids, as is done in Ref. .


text:
# From Mott insulator to ferromagnetic metal: a pressure study of Ca2RuO4

## Abstract

We show that the pressure-temperature phase diagram of the Mott insulator Ca2RuO4 features a metal-insulator transition at 0.5GPa: at 300K from paramagnetic insulator to paramagnetic quasi-two-dimensional metal; at $`T`$ 12K from antiferromagnetic insulator to ferromagnetic, highly anisotropic, three-dimensional metal. We compare the metallic state to that of the structurally related $`p`$-wave superconductor Sr2RuO4, and discuss the importance of structural distortions, which are expected to couple strongly to pressure.

The Mott transition from the metallic to the insulating state, driven by a coulomb interaction, is one of the most dramatic and fundamental many-body effects in condensed matter physics. In recent years there has been a resurgence of interest in the Mott transition due on the one hand to the discovery of materials such as the high-$`T_c`$ cuprates and the manganites in which the metal-insulator transition plays an essential role imada , and on the other to new insights gained from theoretical developments such as the dynamical mean field theory georges  and the recognition of the importance of orbital ordering at the metal-insulator transition imada .

The system Ca2-xSrxRuO4 has recently been investigated in detail and shown to have a metal-insulator transition at 0K at $`x0.2`$ SN\_CSRO . This is a layered perovskite system, isostructural to the cuprate La2-xSrxCuO4, with one end-member, Sr2RuO4, being a highly anisotropic metal with an unconventional (probably $`p`$-wave) superconducting ground state SRO\_YM ; maeno\_pt , while the other, Ca2RuO4, is an antiferromagnetic (AF) insulator for $`T<113`$K, a paramagnetic (PM) insulator for $`113`$K$`<T<360`$K, and a “bad metal” for $`T>360`$K SN\_CSRO ; cao\_1 .

This system is of great interest for the following reasons:

Firstly, the overwhelming majority of work on the metal-insulator transition has been done on $`3d`$-transition metal compounds (mainly oxides) such as V2O3 imada . $`3d`$-orbitals are comparatively compact and close to the ionic core, giving rise naturally to small bandwidths $`W`$, and large Coulomb repulsion $`U`$ between $`3d`$-electrons on the same site. In contrast, $`4d`$ metal-insulator systems are rare, but they should offer a different viewpoint because the more extended $`4d`$-orbitals will tend to have larger effective correlation energy $`U/W`$, but greater susceptibility to orbital ordering via structural distortions around the $`4d`$-site. Indeed the phase diagram of Ca2-xSrxRuO4 shows a series of structural transitions that are intimately related to changes in the electronic state SN\_CRO .

Secondly, although the cuprates show superconductivity in close proximity to a Mott insulating state, a clear connection between the metal-insulator transition and superconductivity has not been made experimentally, and there are many unresolved issues associated with the physics of the cuprates near the Mott insulating state. As a quasi-two-dimensional (Q2D) system bridging the gap between a Mott insulator and an unconventional superconductor, the Ca2-xSrxRuO4 system promises to shed light on this question.

Finally, detailed understanding of how $`p`$-wave superconductivity arises in Sr2RuO4 is lacking. $`p`$-wave superconductivity is assumed to be mediated by magnetic fluctuations at low $`q`$, normally associated with proximity to a ferromagnetic (FM) state at 0K (i.e. to a FM quantum critical point), but Sr2RuO4 does not itself appear to be close to FM order. The doped system, however, at $`x0.5`$ is on the verge of FM order at 0K SN\_CRO , although long-range FM order does not develop for any value of $`x`$.

Pressure is generally to be preferred to chemical doping to tune the properties of materials because pressure does not introduce disorder. Disorder can be extremely hostile to unconventional superconducting states mackenzie\_prl  and it also introduces enormous complications in studying the Mott transition. In this paper we show that application of very modest pressures to pure Ca2RuO4 transforms it from an insulator to a highly anisotropic metal with a FM metallic ground state, and we follow the evolution of this ground state as a function of pressure. Our results shed new light on the nature of ferromagnetism in the layered ruthenates, and open up a new $`4d`$-system for the study of the metal-insulator transition in the absence of disorder.

We used single crystal Ca2RuO4 with an essentially stoichiometric oxygen content (S-CRO) SN\_CRO . Resistivity under pressure $`P`$ was measured by a standard four-probe method, with a ring contact geometry for the $`c`$-axis resistivity $`\rho _c`$. The electrodes were made by a silver evaporation technique. The typical sample size was $`0.15\times 0.5\times 0.04`$ mm3 for in-plane resistivity $`\rho _{ab}`$, and $`0.35\times 0.35\times 0.04`$ mm3 for $`\rho _c`$. Pressures up to 8GPa were generated in a cubic anvil device with an equal volume mixture of Fluorinert FC70 and FC77 (3M Co. St. Paul, MN, USA) as a pressure-transmitting medium. In this device quasi-hydrostaticity can be generated by isotropic movement of six anvil tops even after the fluid medium vitrifies at low temperature and high pressure.

Figure 1 shows $`\rho _{ab}`$ and $`\rho _c`$ at 300K as a function of pressure up to 8 GPa. Both $`\rho _{ab}`$ and $`\rho _c`$ drop discontinuously at 0.5 GPa, indicating an insulator-metal transition. The discontinuity indicates that the transition is of first order and is thus most likely accompanied by a structural change friedt\_unpub .

In the insulating state both $`\rho _{ab}`$ and $`\rho _c`$ decrease isotropically and monotonically with increasing pressure. Following the discontinuous drop at 0.5 GPa however, while $`\rho _{ab}`$ continues to fall, $`\rho _c`$ begins to increase. The anisotropy ratio, $`\rho _c`$/$`\rho _{ab}2`$ in the insulating phase, jumps to about 100 in the metallic phase. Thus, the Q2D nature of the metallic state is comparable with high $`T_c`$ systems or Sr2RuO4.

Figure 2 shows $`\rho _{ab}`$ and $`\rho _c`$ as a function of temperature at several fixed pressures. At ambient pressure, both $`\rho _{ab}`$ and $`\rho _c`$ show an activation-type increase. In contrast, above 0.5 GPa Ca2RuO4 shows a Q2D-metallic nature for temperatures above $``$25K. That is, $`\rho _{ab}(T)`$ shows metallic (d$`\rho _{ab}`$/d$`T>`$0) conduction whereas $`\rho _c`$ has a negative slope indicating nonmetallic conduction. The anisotropy ratio rises to $``$1000 at 25K. The nonmetallic behaviour of $`\rho _c`$ down to 25K in pressurised Ca2RuO4 is in sharp contrast with the metallic $`\rho _c`$ seen below 130K in Sr2RuO4. The rise in $`\rho _c`$ with increasing pressure at 300K also contrasts with simple ideas of metallic conduction, in which the bandwidth would increase with pressure; but the same effect is seen in Sr2RuO4 at $`T>130`$K yoshida , so this seems to be a signature of incoherent $`c`$-axis transport in the ruthenates. We infer that above 25K the $`c`$-axis conduction in metallic Ca2RuO4is mainly described by a tunnelling mechanism as in the high $`T_c`$’s.

As the temperature falls below $``$25K $`\rho _c(T)`$ peaks and then decreases rapidly, indicating a change from two-dimensional to three-dimensional metallic conduction. At the temperature $`T^{}`$ corresponding to the kink in the $`\rho _c(T)`$ curves, $`\rho _{ab}(T)`$ also exhibits an abrupt change in slope, falling rapidly to a residual resistivity of $`3\mu \mathrm{\Omega }`$ cm. This is a very low resistivity for an oxide metal, indicating that our samples are of high quality, but it is still somewhat higher than is required to observe superconductivity in Sr2RuO4 mackenzie\_prl . The behaviour near 25K is reminiscent of a phase transition. In order to clarify this point, the field ($`H`$) and temperature ($`T`$) dependence of the magnetization have thus been measured under pressures of 0.1 MPa, 0.2, 0.5, 0.7, and 0.8 GPa using a piston-cylinder Be-Cu clamp cell with a commercial SQUID magnetometer (Quantum Design, model MPMS) HPcell . We have used several rectangular pieces of sample with total mass 12.5mg for $`H_c`$ and 6mg for $`H_c`$.

Figure 3 shows magnetization ($`M`$) curves at 2 K with $`H_c`$. FM ordering is indicated by the hysteresis in the $`M`$-$`H`$ curves at 0.7 GPa, a pressure at which the system is metallic in the corresponding resistivity measurements. In contrast, no hysteresis is observed in the $`M`$-$`H`$ curves at 0.1 MPa and 0.2 GPa where Ca2RuO4 is in the AF insulating phase. The remnant magnetisation is 0.09 $`\mu _\mathrm{B}`$ at 0.7 GPa, growing to 0.125 $`\mu _\mathrm{B}`$ at 0.8 GPa note2 , while the coercive force of $`H_\mathrm{c}`$ 50mT, which is defined as the field necessary to restore zero magnetisation, increases only weakly.

The volume fraction of the AF insulating phase ought to decrease with pressure above 0.5 GPa, so the remnant magnetisation originates from the metallic phase. Moreover the small moment compared with the saturated moment 2$`\mu _\mathrm{B}`$ of a localized Ru4+ ion implies itinerant magnetism. The shape of the hysteresis curve, which indicates a typical soft ferromagnet, is different from that for the canted AF case seen for example in the AF insulating phase of Ca2RuO4 with excess oxygen (L-CRO) SN\_CRO . We have also observed a similar shape of hysteresis loop with $`H_c`$ at 0.8 GPa where the remnant magnetisation is 0.08 $`\mu _\mathrm{B}`$ and the coercive force of $`H_\mathrm{c}`$ is $`50`$mT. Comparing the hysteresis measured in $`H_c`$, we do not find anisotropy typical of antiferromagnetism. We therefore deduce that Ca2RuO4 in the pressure-induced metallic phase is an itinerant ferromagnet at low $`T`$, although neutron scattering should be done to confirm this.

Figure 4 summarises the pressure dependence of $`T_\mathrm{N}`$ and $`T_\mathrm{C}`$ obtained from our $`M(T)`$ curves, and of $`T^{}`$ obtained from resistivity. $`T_\mathrm{N}`$ is almost independent of pressure up to 0.5GPa, where it abruptly vanishes and is replaced by FM order below $``$12K. The FM $`T_\mathrm{C}`$ then rises with increasing pressure, becoming comparable to $`T^{}`$ near 1GPa. Thus, the remarkable drop in $`\rho (T)`$ at $`T^{}`$ is most probably due to reduced magnetic scattering at a FM ordering transition. $`T^{}`$ continues to rise, reaching 25K at $``$5 GPa, then it decreases gradually. Extrapolation of the $`T^{}`$ vs $`P`$ curve suggests that the FM order may be completely suppressed at $``$15GPa.

From $``$360 to 113K the 0.5 GPa transition is from a PM insulator to a PM-Q2D metal. The PM insulating phase with thermally activated conductivity is important because it shows that Ca2RuO4 is a many-body insulator, as opposed to being a band-insulator with fully opened gaps on the Fermi surface due to AF order. It is also instructive to contrast the PM metallic state with the metallic phase of the manganites, which are always FM and can only be reached by doping. This reflects the difference in on-site repulsion $`U`$ in $`3d`$ and $`4d`$ systems. Mn has a large $`U`$, so the Mn ion in the manganites is in a high-spin configuration, and the strong Hund’s rule coupling requires that the spin of a conduction electron in an $`e_g`$ orbital must be parallel to the localised $`t_{2g}`$ spins both on the ion it is hopping from and the one it is hopping to (so-called “double exchange”). In the ruthenates in contrast lower $`U`$ produces a low-spin configuration, and there is no constraint on the spin-direction of the hopping electron in the metallic state.

It is tempting to interpret the phase diagram of Ca2RuO4 as a classic Mott transition to a metal with a FM ground state Fresard . But this might be too simple. Analogy with Ca2-xSrxRuO4 shows that the physics is not simply driven by changes in $`U/W`$, but also the orbital degeneracy of the Ru4+ $`t_{2g}`$ levels may abruptly change due to structural transitions. These have been well characterised by neutron and x-ray scattering SN\_CSRO ; MB\_CSRO . Pure Sr2RuO4 is tetragonal, but at the Ca rich end the structure is distorted by rotations of the oxygen octahedra about the $`c`$ axis, tilts about an in-plane axis, and large variations in the degree of ‘flattening’ of the octahedra along the $`c`$ axis. In pure Ca2RuO4 the $`T>360`$K metallic phase is L-Pbca, with a smaller volume and weaker tilt and flattening distortions than the $`T<360`$K insulating S-Pbca phase. Pressure should stabilise the low volume phase, so we expect the structural transition associated with the first-order insulator-metal transition to be from L-Pbca to S-Pbca friedt\_unpub . Fang and Terakura Fang  have performed electronic structure calculations that take account of the distortions. They find that flattening, which lowers the $`d_{xy}`$ orbital relative to the $`d_{xz}`$ and $`d_{yz}`$ orbitals, favours an insulating AF ground state; the tilt stabilises AF correlations in the metallic state; while rotations favour ferromagnetism. Thus the insulator-metal transition and AF order may couple separately to pressure, the former caused by suppression of flatting, the latter by reduced tilting.

We now turn to a discussion of the FM metallic state above the insulator-metal transition. According to band structure calculations FM order in layered ruthenates is a reflection of rotational distortions Fang ; Hotta , which develop FM coupling by enlarging the density of states at the Fermi energy via band narrowing. This idea is supported by the large FM correlations in Ca2-xSrxRuO4 with $`x=`$ 0.5, in which the octahedra are rotated but not tilted unpub . More generally mazin , ruthenates show a tendency to ferromagnetism due to an unusually large admixture at the Fermi surface of the oxygen $`p`$-orbitals, which can gain exchange energy from FM order. It will be very interesting to have electronic structure calculations on metallic Ca2RuO4, once the crystal structure is known accurately in the high-pressure phase. It will also be interesting to compare the properties of the itinerant-FM Q2D-metal with theoretical predictions Hatatani .

Finally, the large range of temperatures for which pressurised Ca2RuO4 exhibits a non-metallic $`\rho _c`$ is of interest. In Ca2-xSrxRuO4, for 0.5 $`<x<`$ 1.5 where the octahedra are rotated but not tilted, nonmetallic $`\rho _c(T)`$ is also seen over a wide temperature range, while Sr2RuO4, having non-rotated octahedra, exhibits metallic $`\rho _c`$ below 130K. We therefore suspect that the rotational distortion affects not only the magnetic properties but also the two-dimensionality in the electronic transport in the metallic phase.

In conclusion, high-purity Ca2RuO4 is a $`4d`$-system that at 0.5 GPa has a transition from a Mott insulating state to a metal with a FM ground state. The $`P`$-$`T`$ phase diagram of this system is unique and quite different from that of the doped system SN\_CSRO ; HF\_LCRO . Furthermore, our results suggest the existence of a FM quantum critical point at pressures above 10 GPa. We thus have strong interest in measurements at higher pressures as a crucial test of the connection between the $`p`$-wave superconductivity of Sr2RuO4, and the ferromagnetism of metallic Ca2RuO4.

We thank T.Moriya, G.G. Lonzarich, A.P. Mackenzie and D.J. Singh for informative discussions. A part of this work was supported by a Grant-in-Aid for Scientific Research on Priority Areas, from the MECSST of Japan, and by the EPSRC of Great Britain.