tcb@breakable

# Domain Decomposition in space-time for 4D-Var Data Assimilation problem:
a case study on the Regional Ocean Modeling System

R. Cacciapuotia, L. D’Amorea and A. M. Mooreb

(a University of Naples, Federico II, Naples (IT).
b Ocean Sciences Department Institute of Marine Sciences, UC Santa Cruz, (CA-
USA).)

##  1 ROMS

Regional Ocean Modeling System (ROMS) is an open-source, mature numerical
framework used by both the scientific and operational communities to study
ocean dynamics over 3D spatial domain and time interval.
ROMS supports different 4D-Var data assimilation (DA) methodologies such as
incremental strong constraint 4D-Var (IS4D-Var) and dual formulation 4D-Var
using restricted B-preconditioned Lanczos formulation of the conjugate
gradient method (RBL4DVAR) [8]. IS4D-Var and RBL4DVAR search best circulation
estimate in space spanned by control vector and observations, respectively.
IS4D-Var and RBL4DVAR algorithm consist of two nested loop, the outer-loop
involves the module (1), namely nonlinear ROMS (NLROM) solving ROMS equations,
the inner-loop involves modules (2)-(3), namely tangent linear approximation
of ROMS (TLROMS) and adjoint model of ROMS (ADROMS); TLROMS and ADROMS are
used for minimizing 4D-Var functional [3] (see Figures 1-2).
NLROMS is a three-dimensional, free-surface, terrain-following ocean model
that solves the Reynolds-averaged Navier-Stokes equations using the
hydrostatic vertical momentum balance and Boussinesq approximation.
NLROMS computes

$x^{ROMS}(t_{l})=M_{l-1,l}(x(t_{l-1}),f(t_{l}),b(t_{l}))$ (1)

with the state-vector $x^{ROMS}(t_{l})=(T,S,\varsigma,u,v)^{T}$, temperature
$T$, salinity $S$, $(x,y)$ components of vector velocity $u,v$, sea surface
displacement $\varsigma$. $M_{l-1,l}$ represents nonlinear ROMS acting on
$x^{ROMS}(t_{l-1})$, and subject to forcing $f(t_{l})$, and boundary
conditions $b(t_{l})$ during the time interval $[t_{l-1},t_{l}]$.

Minimization of the 4D-Var functional:

$\mathbf{J}^{ROMS}(\delta z)=\frac{1}{2}\delta z\mathbf{B}^{-1}\delta
z+\frac{1}{2}(G\delta z-\mathbf{d})^{T}\mathbf{R}^{-1}(G\delta z-\mathbf{d})$
(2)

where $\delta z$ are the control variable increments, $\mathbf{d}$ is vector
of innovations, $G=(...,H_{l}^{T},...)^{T}$, where $H_{l}$ is the observation
matrix; $\mathbf{R}$ is observation error covariance matrix and $\mathbf{B}$
is covariance matrix of model error, is computed in the inner-loop in Figure
1.
Analysis increment, $\delta z^{a}$, that minimizes 4D-Var function in (2)
corresponds to the solution of the equation $\
\partial\mathbf{J}^{ROMS}/\partial\delta z=0$, and is given by:

$\delta
z^{a}=(\mathbf{B}^{-1}+G^{T}\mathbf{R}^{-1}G)^{-1}G^{T}\mathbf{R}^{-1}\mathbf{d}$
(3)

or, equivalently

$\delta z^{a}=\mathbf{B}G^{T}(G\mathbf{B}G^{T}+\mathbf{R})^{-1}\mathbf{d}.$
(4)

Equation (3) is referred to as the dual form (RBL4DVAR), while (4) is referred
to as the primal form (IS4DVAR). In particular, we define

$\mathbf{K}=\mathbf{B}G^{T}(G\mathbf{B}G^{T}+\mathbf{R})^{-1}$ (5)

as Kalman gain matrix.

TLROMS computes

$\delta x^{ROMS}(t_{l})\simeq M_{l-1,l}u(t_{l-1})$ (6)

where $\delta x^{ROMS}(t_{l})=x^{ROMS}(t_{l})-x^{b}(t_{l})$, $\delta
f(t_{l})=f(t_{l})-f^{b}(t_{l})$, $\delta b(t_{l})=b(t_{l})-b^{b}(t_{l})$, and
$x^{b}(t_{l})$, $f^{b}(t_{)}$, $b^{b}(t_{l})$ are the background of the
circulation, surface forcing and open boundary conditions respectively, and

$u(t_{l-1})=((\delta x^{ROMS})^{T}(t_{l-1}),\delta f^{T}(t_{l}),\delta
b^{T}(t_{l}))^{T}.$

Equation (6) is obtained from first-order Taylor expansion of NLROMS in (1).

ADROMS computes

$u^{*}(t_{l-1})=M_{l-1,l}^{T}p(t_{l})$ (7)

where $u^{*}(t_{l-1})=(p^{T}(t_{l-1}),\delta f^{T}(t_{l}),\delta
b^{*T}(t_{l}))^{T}$ where $p$ is the adjoint state-vector, $\delta f^{T}$ and
$\delta b^{*T}$ are the adjoint of the surface forcing and the open boundary
condition increments.

##  2 DD-4DVarDA in ROMS model

Domain Decomposition (DD) method proposed in [1] is made up of decomposition
of the domain $\Omega\times\Delta$ into subdomains where $\Omega$ is the 3D
spatial domain and $\Delta$ is the time interval, solution of reduced forecast
model and minimization of local 4D-Var functionals (we call it DD-4DVarDA
method).
Relying on the existing software implementation, in the next we describe main
components of DD-4DVarDA method, highlighting the topics that we will address
both on the mathematical problem underlying ROMS and the code implementation
(see steps 1-4). We focus on IS4DVAR formulation described in Section 1.

###  2.1 Decomposition in space and time of the ocean model

Figure 1: A flow chart illustrating IS4D-Var algorithm. Figure 2: The flow
chart illustrating RBL4DVAR algorithm.

Decomposition of spatial domain $\Omega$.
We will consider a 2D decomposition of $\Omega\subseteq\mathbb{R}^{3}$ in x-
and y-direction and denote $\Omega_{xy}$ the spatial domain to decompose.

ROMS uses a parallelization approach that partitions domain $\Omega_{xy}$ into
tiles (see Figure 3) [2]

$\Omega_{xy}=\bigcup_{i=0}^{N_{sub}-1}tile_{i}$ (8)

where $N_{sub}=NtileI\times NtileJ$; $NtileI$ and $NtileJ$ are the number of
tiles set in the input file in x- and y-direction, respectively.
We denote by $HI$ and $HJ$ the overlapping tiles regions (i.e ghost or halo
area in ROMS111In ROMS, the halo area would be two grids points wide unless
the MPDATA advection scheme is used, in which case it needs three.) in x- and
y-direction, respectively.

Step 0: DD of $\Omega$ in ROMS.
In our study we will assume the decomposition available in ROMS, as given in
(8).

Decomposition of time interval $\Delta$.
ROMS does not yet implement decomposition in time direction.

Step 1: DD of $\Delta$ in ROMS.
In our study we aim introducing a decomposition of time interval $\Delta$ into
$N_{t}$ intervals:
$\Delta=\bigcup_{k=1}^{N_{t}}\Delta_{k}:=\bigcup_{k=1}^{N_{t}}[t_{\bar{s}_{k-1}},t_{\bar{s}_{k-1}+N_{k}}],$
(9) where $N_{k}=|D(\Delta_{k})|$ are respectively the number of subdomains of
$[0,T]$ and of time $t_{l}\in\Delta_{k}$ such that
$\sum_{k=1}^{N_{t}}N_{k}-(N_{t}-1)=N$,
$\bar{s}_{k-1}:=\sum_{j=1}^{k-1}N_{j}-(k-1)$ and $\bar{s}_{0}:=0$.

Figure 3: A tiled grid in xy-plane with some internal ROMS parameters.

Ocean model reduction. ROMS allows each tile (or subdomain, see Figure 3) to
compute local solutions of TLROMS and ADROMS.
For $i=0,1,\ldots,N_{sub}-1$ and $k=1,\ldots,N_{t}$ local222Let
$x\in\mathbb{R}^{N_{p}}$ and $y\in\mathbb{R}^{N_{p}\times N}$ be vectors, for
simplicity of notations, we refer to $x_{i}$ as a restriction of $x$ to
$\Omega_{i}$, i.e. $x_{i}\equiv x/\Omega_{i}$ and $x_{i,k}\equiv
x/(\Omega_{i}\times\Delta_{k})$, similarly for matrix
$A\in\mathbb{R}^{N_{p}\times N_{p}}$, i.e. $A_{i}\equiv A/\Omega_{i}$ and
$A_{i,k}\equiv A/(\Omega_{i}\times\Delta_{k})$, according the description in
[6]. TLROMS on local domain $tile_{i}$ computes

$\delta x_{i}^{ROMS}(t_{l})\simeq M_{i,(l-1,l)}u_{i}(t_{l-1})$ (10)

and local ADROMS on local domain $tile_{i}$ computes

$u_{i}^{*}(t_{l-1})=M_{i,(l-1,l)}^{T}p_{i}(t_{l}),$ (11)

where $x_{i}$, $u_{i}$, $u_{i}^{*}$, $M_{i}$, $p_{i}$ are the restriction on
$tile_{i}$ of variables $x$, $u$, $u^{*}$, $M$ and $p$ and $M_{i,(l-1,l)}$ is
discrete model from $t_{l-1}$ to $t_{l}$.
DD-4DVar method introduces model reduction by using the background $x^{b}$ as
local initial values. For $n=0,1,\ldots,\bar{n}$ (outer loop of DD–4DVAR
method [1]) do: for $k=1,\ldots,N_{t}$, posed $x_{i,k}^{0}\equiv x_{i,k}^{b}$
$\forall i=0,1,\ldots,N_{sub}-1$, we let $x_{i,k}^{{M}_{i,k}}$ be the solution
of the local model

$(P_{i,k}^{{M}_{i,k},n})_{i=0,1\ldots,N_{sub}-1,r=1,\ldots,N_{t}}\
:\quad\left\\{\begin{array}[]{ll}x_{i,k}^{\mathcal{M}_{i,k},n}=M_{i,k}x_{i,k-1}^{n}+b_{i,k},\\\
x_{i,k-1}^{n}=x_{i,k}^{\mathcal{M}_{i,k},n},\\\
x_{i,k-1}^{n}/HI=x_{i_{I},k-1}^{n}/HI,\ \ (\ref{co}.1)\\\
x_{i,k-1}^{n}/HJ=x_{i_{J},k-1}^{n}/HJ,\ (\ref{co}.2)\end{array}\right.$ (12)

where $i_{I}=0,\ldots,n_{I}-1$, $i_{J}=0,\ldots,n_{J}-1$, $n_{I}$ and $n_{J}$
are respectively numbers of adjacent tiles in x- and y-direction,
$b_{i,r}^{k}$ and $M_{i}^{r}$ are respectively background on
$tile_{i}\times\Delta_{r}$, the vector accounting boundary conditions of
$tile_{i}$ and the restriction in $tile_{i}$ of the matrix in (1) that is

$M_{k}\equiv
M_{\bar{s}_{r-1},\bar{s}_{r}}:=M_{\bar{s}_{r-1},\bar{s}_{r-1}+1}\cdots
M_{\bar{s}_{r}-1,\bar{s}_{r}}.$ (13)

In the following, we neglect the dependency on outer loop iteration $n$ of
DD–4DVAR method. We underline that local TLROMS and ADROMS in (6) and (7) are
obtained by using MPI exchange for boundary conditions, regardless of local
solution on overlap area, namely they do not consider overlapping tiles
conditions in (12.1) and (12.2). Consequently, we need to modify local TLROMS
and ADROMS in (10) and (11) taking into account of (12.1) and (12.1)
conditions. More precisely,

Step 2: TLROMS.
for $k=1,\ldots,N_{t}$, $\forall i=0,1,\ldots,N_{sub}-1$, local TLROMS on
$tile_{i}\times\Delta_{k}$ will be modified such that $\delta x_{i,k}\simeq
M_{i,k}u_{i,k-1}+\theta_{I}(u_{i,k-1})+\theta_{J}(u_{i,k-1})$ (14) where
$\begin{array}[]{cc}\theta_{J}(u_{i,k-1}):=&\sum_{i_{J}=1}^{n_{J}}\gamma_{i_{J}}(M_{i,k}/HJ\cdot
u_{i,k-1}/HJ-\\\ &M_{k}/HJ\cdot u_{i,k-1}/HJ)\end{array}$ (15) and
$\begin{array}[]{cc}\theta_{I}(u_{i,r-1}):=&\sum_{i_{I}=1}^{n_{I}}\gamma_{i_{I}}(M_{i}^{r}/HI\cdot
u_{i,r-1}/HI-\\\ &M_{i}^{r}/HI\cdot u_{i,r-1}/HI)\end{array}$ (16) are
overlapping vectors in x- and y-direction, respectively; where parameters
$\gamma_{i_{I}}$, $\gamma_{i_{J}}$ denote weights.

Similarly, we need to modify local ADROMS in (11). More precisely,

Step 3: ADROMS.
for $k=1,\ldots,N_{t}$, $\forall i=0,1,\ldots,N_{sub}-1$, local ADROMS in (11)
on local domain $tile_{i}\times\Delta_{k}$ will be modified such that
$u_{i,k-1}=(M_{i,r})^{T}\cdot
p_{i,k}+\theta_{J}(p_{i,k-1})+\theta_{I}(p_{i,k-1})$ (17) where $\theta_{J}$
and $\theta_{I}$ are defined in (15) and (16).

Figure 4: A region within the domain

4D-VAR DA Operator Reduction.
Operator reduction involves the inner loop in Figure 1.
ROMS [2] computes and minimizes the operator

$\mathbf{J}_{tile_{i}}^{ROMS}:=\mathbf{J}^{ROMS}/tile_{i}$ (18)

where $\mathbf{J}^{ROMS}$ is defined in (2) and $\mathbf{J}^{ROMS}/(tile_{i})$
is the IS4D-Var functional in each tile $tile_{i}$ with MPI exchange of
boundary conditions, $\forall i=0,1,...,N_{sub}-1$.

Local 4D–VAR DA functional in [1] is defined as follows

$\mathbf{J}_{tile_{i}\times\Delta_{r}}^{DD-4DVar}:=\mathbf{J}_{i,k}(x_{i,k})=\mathbf{J}(x_{i,k})/(tile_{i}\times\Delta_{k})+\mathcal{O}_{IJ}(x_{i,k})$
(19)

where

$\mathbf{J}(x)=\alpha\|x-x^{b}\|_{\textbf{B}^{-1}}^{2}+\|Gx-y\|_{\textbf{R}^{-1}}^{2},$
(20)

is the DD-4DVar functional in [1],

$\mathcal{O}_{IJ}(x_{i,k})=\sum_{i_{I}=1}^{n_{I}}\beta_{i_{I},k}\cdot\|x_{i,k}/HI-
x_{i_{I},k}/HI\|_{\textbf{B}_{i_{I}}^{-1}}^{2}+\sum_{i_{J}=1}^{n_{J}}\beta_{i_{J},k}\cdot\|x_{i,k}/HJ-
x_{i_{J},k}/HJ\|_{\textbf{B}_{i_{j}}^{-1}}^{2}$ (21)

is the overlapping operator on overlapping tiles region $HI$ and $HJ$, and

$\mathbf{J}_{i,k}(x_{i,k})/(tile_{i}\times\Delta_{k})=\alpha_{i,k}\cdot\|x_{i,k}-x_{i,k}^{{M}_{i,k}}\|_{\textbf{B}_{i}^{-1}}+\|{G}_{i,r}x_{i,k}-y_{i,k}\|_{\textbf{R}_{i}^{-1}}^{2}$
(22)

is the restriction of $\mathbf{J}$ on $tile_{i}\times\Delta_{k}$, where
$x^{b}$ is background, $y$ is observations vector in $\Delta$, $y_{i,k}$ is
observations vector in $\Delta_{k}$; $\mathbf{B}_{i}$, $\mathbf{B}_{i_{I}}$,
$\mathbf{B}_{i_{J}}$ are respectively the restrictions of covariance matrix
$\mathbf{B}$ to $tile_{i}$, $HI$ and $HJ$; $G_{i,r}$, $\mathbf{R}_{i}$ are
restriction of matrices $G$ and $\mathbf{R}$ to $tile_{i}$ in $\Delta_{k}$.
Parameters $\alpha_{i,r}$, $\beta_{i_{I}}$ and $\beta_{i_{J}}$ in (21) denotes
regularization parameters. We let
$\alpha_{i,k}=\beta_{i_{I},k}=\beta_{i_{J},k}=1$, $\forall
i_{I}=1,\ldots,n_{I}$ and $\forall i_{J}=1,\ldots,n_{J}$.
Incremental formulation of local 4D–VAR DA functional in (19) is

${J}_{i,k}(\delta z_{i,k})={J}(\delta
z_{k})/(tile_{i}\times\Delta_{k})+\mathcal{O}_{IJ}(\delta z_{i,k})$

where $\delta z_{i,k}$ are the control variable increments in
$tile_{i}\times\Delta_{k}$.

Step 4: IS4D-Var.
From (20) local IS4D-Var function in (18) can be written
$\mathbf{J}_{loc}^{ROMS}=\mathbf{J}_{tile_{i}\times\Delta}^{DD-4DVar}-\mathcal{O}_{IJ}$
(23) Consequently, we need to add overlapping operator in order to enforce the
matching of local solutions on the overlapping tiles in each time interval.

Figure 5:  Placement of variables on an Arakawa C grid.  Figure 6:  ROMS
directory structure.

###  2.2 DD-4DVarDA in ROMS code

In Figure 6 the ROMS directory structure is shown. We focus on ROMS folder, in
particular, on its folders: Tangent and Adjoint.

  1. 1.

ROMS.
We need to modify routines in ROMS folder implementing decomposition as in (9)
(see step 1).
Decomposition of time interval involves modification of initial conditions of
nonlinear model, in TLROMS and ADROMS.
Routines involving initialization are (see Figure 6):

    1. (a)

initial routine in Nonlinear folder: initializes all model variables (it is
called in main3d).

Note 1.1 (see step 1)
Initial routine initializes all model variable before calling main3d routine
(in Nonlinear folder). Main3d routine is the main driver solving nonlinear
ROMS model (background). DD in x- and y- directions is not applied for
computing background, i.e. there are not MPI communications among processes.
This means that background is computed without using DD in space, consequently
we not apply DD in time. Moreover, some values of $u,v$ (components velocity)
are not set to zero at the end of time step, because its values on some grid
points are necessary for the next time step.

    2. (b)

tl_initial routine in Tangent folder: initializes all tangent model variables
(it is called in i4dvar).

Note 1.2 (see step 1).
We consider tl_initial routine. Parts involving initializations are: • line
123: Initializes time stepping indices and counter. • line 162: initialization
time.

    3. (c)

ad_initial routine in Adjoint folder: initializes all adjoint model
variables(it is called in i4dvar).

Note 1.3 (see step 1).
We consider ad_initial routine. Parts involving initializations are: • line
113: Initializes time stepping indices and counter. • line 152: initialization
time.

Main actions to apply DD in time in ROMS are described below.

Action 1.1 (see step 1).
(a) Modify initials routines by adding MPI communications for initial
conditions to tangent and adjoint routines in each time interval.
(b) Decompose time interval (RunInterval variable in routines). One possible
way forward: • We need to use OMP threadprivate directive for replicating
variables related to time interval such that each thread has its own copy.
(see mod_parallel routine in Module folder at line 51 related to DD in space).
Allocate_routine in mod_parallel allocates tiles; equally we can allocate
local time interval and related variables.
• We can add allocation and OMP threadprivate directive of local time
intervals in mod_parallel routine. We can define first time interval
(first_time_interval) and last time interval (last_time_interval) and add a
for loop, after for loop involving tiles, started from first time interval up
to last time interval adding a time step (dt) defined in driver.
– tl_main3d: at line 277 stars for loop involves tiles. – ad_main3d: at line
629 stars for loop involves tiles. • Identify tangent and adjoint variables
for MPI communications in time. (c) Introduce MPI communications in time. • We
need to split ROMS MPI communicator (OCN_COMM_WORLD) for obtaining MPI
communicators needed to communications among processes related to same spatial
sub domain but different time intervals.
By splitting ROMS MPI communicator we obtained a new MPI communicator namely
one single communicator for the current process.

Action 1.2 (see step 1).
Step 1 involves: • tl_main3d routine: at line 142 starts while loop on time
interval (RunInterval) by increasing the step time (my_steptime). Moreover, at
lines 278 and 280 it calls tl_set_massflux and tl_rho_eos, i.e. the routines
we need to modify (see Action 2.1,2,3,4). Consequently, we probably need to
introduce in tl_main3d MPI communications in time after the tl_set_massflux
and tl_rho_eos routines. • ad_main3d routine: at line 177 starts while loop on
time interval (RunInterval) by increasing the step time (my_steptime).
Moreover, at lines 632 and 634 it calls ad_rho_eos and ad_set_massflux, i.e.
the routines we need to modify (see Action 3.1 and 3.2). Consequently, we
probably need to introduce in ad_main3d MPI communications in time after the
ad_set_massflux and ad_rho_eos routines.

  2. 2.

TLROMS.
We need to modify TLROMS as in (14) (see step 2), namely we need to compute
the overlapping vector in (15) and (16) and add them to tangent variables.

Action 2.1 (see step 2).
We could consider another inner loop (inside the loop over $m$ in Figure 1)
over index $n$ and initial approximation of solution on local tiles at $n=0$.
For each iteration we need local solution on tile adjacent to each tile, this
means using MPI exchange of information between adjacent tiles (in two
direction see Figure 5).

We note that primitive equations of motion [5] are written in flux form
transformed using orthogonal curvilinear coordinates ($\xi,\eta$) (see Figure
5).
We consider tl_main3d routine. In tangent folder (see Figure 6) this routine
is the main driver of TLROMS configurated as a full 3D baroclinic ocean model.
tl_main3d routine calls the following subrountines.

    1. (a)

tl_set_massflux calls tl_set_massflux_tile.
tl_set_massflux_tile: computes “in situ” tangent linear horizontal mass flux.

Action 2.2 (see step 2).
Taking into account (15) and (16) we need to modify the code starting from
line 155.

    2. (b)

tl_rho_eos calls tl_rho_eos_tile.
tl_rho_eos_tile: computes “in situ” the density and other quantities
(temperature, salinity,…).

Action 2.3 (see step 2).
Parts to be modified: i. Line 505: computes ”in situ” density anomaly; ii.
Line 591: computes ”in situ” Brunt-Vaisala frequency; iii. Line 1158: computes
”in situ” Brunt-Vaisala frequency;

    3. (c)

tl_omega: computes vertical velocity (no modifications because DD is only
horizontal, there are not overlap region along vertical direction).

  3. 3.

ADROMS.
Similarly to TLROMS, we need to modify ADROMS as in (17) (see step 3), namely
we need to compute and add the overlapping vector in (15) and (16) to adjoint
variables.
We consider ad_main3d routine. In Adjoint folder (see Figure 6) this routine
is the main driver of ADROMS configurated as a full 3D baroclinic ocean model.
ad_main3d routine calls the following subrountines.

    1. (a)

ad_rho_eos: computes “in situ” density and other associated quantities.

Action 3.1 (see step 3).
Parts that should be modified: i. Lines 797 and 1758: compute ”in situ”
adjoint Brut-Vaisala frequency at horizontal points. ii. Line 1070: computes
”in situ” adjoint density anomaly.

    2. (b)

ad_set_mass_flux: compute “in situ” adjoint horizontal mass fluxes.

Action 3.2 (see step 3).
Part that should be modified: i. Line 201: computes ”in situ” adjoint
horizontal mass fluxes

    3. (c)

ad_set_avg: accumulates and computes output time-averaged adjoint fields.
(probably no modifications are needed).

  4. 4.

IS4D-Var.
We need to modify routines in ROMS folder as in (23) (see step 4).

Action 4.1 (see step 4).
We need to modify IS4DVAR cost function in (2) to take in account overlap
region i.e. halo region.

##  3 ROMS test cases

We have installed ROMS on the high-performance hybrid computing architecture
of the Sistema Cooperativo Per Elaborazioni scientifiche multidiscipliari data
center, located in the University of Naples Federico II. More precisely, the
HPC architecture is made of 8 nodes, consisting of distributed memory DELL
M600 blades connected by a 10 Gigabit Ethernet technology. Each blade consists
of 2 Intel<EMAIL_ADDRESS>GHz quadcore processors sharing 16 GB RAM memory for a
total of 8 cores/blade and of 64 cores, in total.
We consider test cases referred to IS4D-Var and RBL4DVAR configured for the
U.S. west coast and the California Current System (CCS). This configuration,
referred to as WC13, has 30 km horizontal resolution, and 30 levels in the
vertical and data assimilated every 7 days during the period July 2002 – Dec.
2004.

  * •

IS4DVAR test cases.

To run this application I need to take the following steps.

    * –

I customize build_roms.csh script in IS4DVAR folder and define the path to the
directories where all project’s files are kept.

      * *

I set USE_MY_LIBS to yes, then I modify source code root file
Complilers/my_build_paths.sh and edit the appropriate paths for the desired
compiler.

    * –

I customize the ROMS input file (roms_wc13_2hours.in/roms_wc13_daily.in) and
specify the appropriate values of tile numbers in the I-direction and
J-direction i.e. NtileI=2 and NtileJ=4 (see Figure 3).

    * –

I modify ARmake.inc in Lib folder to create the libraries for serial and
parallel ARPACK.

    * –

I customize Linux-ifort.mk (Fortran compiler is ifort) and define library
locations.

    * –

I create PBS scripts to run the test cases.

Figures 8 and 7 show results obtained applying DD in space by considering
$NtileI\times NtileJ$ tiles and input files roms_wc13_2hours.in and
roms_wc13_daily.in, respectively. They show total cost function $J$ (black
curve), observation cost function $Jo$ (blue curve), and background cost
function $Jb$ (red curve) and the theoretical minimum value
$J_{min}=n_{obs}/2$ (dashed black line) plotted on a $log_{10}$ scale. We
observe that the solution has more or less converged by 25 iterations.
We will refer to the following quantities: Nouter number of outer loops,
Ninner inner loop, $T_{Casei}^{node}$ elapsed time for each node in seconds,
$T_{max}$ and $T_{min}$ maximum and minimum elapsed time among nodes in
seconds.
We consider:

    * –

case 1: Nouter=1 and Ninner=25;

    * –

case 2: Nouter=1 and Ninner=50;

    * –

case 3: Nouter=2 and Ninner=25;

    * –

case 4: Nouter=2 and Ninner=50;

In Tables 1, 2 and 3, 4 we report elapsed time for each process considering
roms_wc13_2hours and roms_wc13_daily.in input file, respectively.

$Node$ | $T_{Case1}^{node}$ | $T_{Case2}^{node}$ | $T_{Case3}^{node}$ | $T_{Case4}^{node}$
---|---|---|---|---
$0$ | $1.56\times 10^{3}$ | $3.24\times 10^{3}$ | $3.124\times 10^{3}$ | $6.42\times 10^{3}$
$1$ | $1.58\times 10^{3}$ | $3.29\times 10^{3}$ | $3.17\times 10^{3}$ | $6.49\times 10^{3}$
$2$ | $1.58\times 10^{3}$ | $3.29\times 10^{3}$ | $3.17\times 10^{3}$ | $6.50\times 10^{3}$
$3$ | $1.58\times 10^{3}$ | $3.28\times 10^{3}$ | $3.16\times 10^{3}$ | $6.50\times 10^{3}$
$4$ | $1.58\times 10^{3}$ | $3.29\times 10^{3}$ | $3.17\times 10^{3}$ | $6.50\times 10^{3}$
$5$ | $1.58\times 10^{3}$ | $3.29\times 10^{3}$ | $3.16\times 10^{3}$ | $6.50\times 10^{3}$
$6$ | $1.58\times 10^{3}$ | $3.29\times 10^{3}$ | $3.16\times 10^{3}$ | $6.50\times 10^{3}$
$7$ | $1.58\times 10^{3}$ | $3.28\times 10^{3}$ | $3.17\times 10^{3}$ | $6.49\times 10^{3}$
Table 1: Elapsed time in seconds for each process related to roms_wc13_2hours input file. | $T_{max}$ | $T_{min}$
---|---|---
${Case1}$ | $1.58\times 10^{3}$ | $1.56\times 10^{3}$
$Case2$ | $3.29\times 10^{3}$ | $3.24\times 10^{3}$
${Case3}$ | $3.17\times 10^{3}$ | $3.12\times 10^{3}$
${Case4}$ | $6.50\times 10^{3}$ | $6.42\times 10^{3}$
Table 2: Maximum and minimum elapsed time in seconds among nodes related to roms_wc13_2hours.in input file. $Node$ | $T_{Case1}^{node}$ | $T_{Case2}^{node}$ | $T_{Case3}^{node}$ | $T_{Case4}^{node}$
---|---|---|---|---
$0$ | $1.41\times 10^{3}$ | $2.95\times 10^{3}$ | $2.79\times 10^{3}$ | $5.97\times 10^{3}$
$1$ | $1.41\times 10^{3}$ | $2.96\times 10^{3}$ | $2.83\times 10^{3}$ | $5.99\times 10^{3}$
$2$ | $1.41\times 10^{3}$ | $2.97\times 10^{3}$ | $2.82\times 10^{3}$ | $5.99\times 10^{3}$
$3$ | $1.41\times 10^{3}$ | $2.97\times 10^{3}$ | $2.82\times 10^{3}$ | $6.00\times 10^{3}$
$4$ | $1.41\times 10^{3}$ | $2.97\times 10^{3}$ | $2.82\times 10^{3}$ | $5.98\times 10^{3}$
$5$ | $1.41\times 10^{3}$ | $2.96\times 10^{3}$ | $2.82\times 10^{3}$ | $5.99\times 10^{3}$
$6$ | $1.41\times 10^{3}$ | $2.96\times 10^{3}$ | $2.82\times 10^{3}$ | $6.00\times 10^{3}$
$7$ | $1.41\times 10^{3}$ | $3.00\times 10^{3}$ | $2.82\times 10^{3}$ | $6.00\times 10^{3}$
Table 3: Elapsed time in seconds for each process releted to roms_wc13_daily.in input file. | $T_{max}$ | $T_{min}$
---|---|---
${Case1}$ | $1.41\times 10^{3}$ | $1.41\times 10^{3}$
$Case2$ | $3.00\times 10^{3}$ | $2.95\times 10^{3}$
${Case3}$ | $2.83\times 10^{3}$ | $2.79\times 10^{3}$
${Case4}$ | $5.97\times 10^{3}$ | $5.97\times 10^{3}$
Table 4:  Maximum and minimum elapsed time in seconds among nodes related to
roms_wc13_daily.in input file.

  * •

RBL4DVAR.

To run this application I need to take the following steps.

    * –

I customize build_roms.csh scripts in RBL4DVAR folders and define the path to
the directories where all project’s files are kept.

    * –

I customize the ROMS input file (roms_wc13_2hours.in/roms_wc13_daily.in) and
specify the appropriate values of tile numbers in the I-direction and
J-direction i.e. NtileI=2 and NtileJ=4 (see Figure 3).

    * –

I create a PBS script to run the test cases.

  * •

IS4D-Var vs RBL4DVAR.

IS4D-Var and RBL4DVAR algorithms are based on a search for the best
circulation estimate in the space spanned by the model control vector and in
the dual space spanned by the observations, respectively. We note that
${n_{obs}}<<N_{p}$, where ${n_{obs}}$ is number of observations and $N_{p}$ is
dimension of model control vector. Hence, the dimension of observation space
is significantly smaller than the model control space, the dual formulation
can reduce both memory usage and computational cost. Consequently, it appears
that the dual formulation should be an easier problem to solve because of the
considerably smaller dimension of the space involved, but it has practical
barriers to convergence.
Consequently, we consider three applications of dual formulation:

    1. 1.

using the standard $R^{-1/2}$ preconditioning with a conjugate gradient
(RBL4DVAR) method;

    2. 2.

using the $R^{-1/2}$ preconditioning and minimum residual algorithm (MINRES);

    3. 3.

using a restricted B-preconditioned conjugate gradient (RPCG) approach.

Figure 9 and 10 shows the convergence of minimization algorithms RBL4DVAR and
RBL4DVAR, MINRES, and RPCG compared to primal formulation IS4DVAR by
considering $Nouter=1$, $Ninner=50$ and $Nouter=1$, $Ninner=25$. We observe
that the performance of RPCG is superior to both MINRES and RBL4DVAR, in
particular, PCG ensures that RPCG converges at same rate as IS4D-Var.

  * •

Observation impact and observation sensitivity.
The observations assimilated into the model are:

    * –

SST satellite: Sea Surface Temperature;

    * –

SSH satellite: Sea Surface Heights;

    * –

Argo floats: hydrographic observations of temperature and salinity.

We consider the time average transport across 37N over the upper 500 m,
denoted by $I_{37N}$, and given by

$I_{37N}(x)=\frac{1}{N}\sum_{i=1}^{N}h^{T}x_{i}$ (24)

where $h$ is a vector with non-zero elements corresponding to the velocity
grid points that contribute to the transport normal to the 37N section shown
in Figure 12, $N$ is the number of time steps during the assimilation
interval, $x_{i}$ is the model state-vector at time $i\Delta t$, and $\Delta
t$ is the model time step. We consider scalar functions of ocean state vector
$I=I(x)$, prior $I^{b}=I(x^{b})$ and posterior ocean state vector
$I^{a}=I(x^{a})$, where $x^{b}$ is background circulation estimate and $x^{a}$
is 4D-Var analysis. The circulation analysis increment $\delta x^{a}(t)$ at
instant time $t$ in interval time $[t_{0},t_{0}+7]$ of each 4D-Var cycle is
given

$\delta x^{ROMS}(t)=x^{a}(t)-x^{b}(t),$ (25)

consequently, transport increment $\Delta I$ during each assimilation cycle
can be expressed as:

$\Delta I=I^{a}-I^{b}=I(x^{b}+\delta x)-I^{b}\simeq\delta x^{T}(\partial
J/\partial x),$ (26)

using the tangent linear assumption in (6), (5) and (4), the circulation
analysis increment is given by

$\delta x^{ROMS}(t)=M(t,t_{0})\tilde{\mathbf{K}}d,$ (27)

where $M(t,t_{0})$ represents the perturbation tangent linear model for the
time interval $[t_{0},t]$, $\tilde{\mathbf{K}}$ is an approximation of
$\mathbf{K}$ defined in (5) and $\mathbf{d}$ is innovation vector. Then,
increment defined in (26) become

$\Delta I\simeq d^{T}\tilde{K}^{T}M^{T}(\partial I/\partial x),$ (28)

where $\tilde{K}$ represents adjoint of Kalman gain matrix.
Consequently, 37N time increment of averaged transport in (24) is:

$\Delta
I_{37N}\simeq\frac{1}{N}d^{T}\tilde{K}^{T}\sum_{l=1}^{N}(M_{l})^{T}h=d^{T}g=d^{T}(g^{x}+g^{f}+g^{b})$
(29)

or

$\Delta I_{37N}\simeq\sum_{l=1}^{n_{obs}}(y_{l}-H_{l}(x^{b}(t)))g_{l},$ (30)

namely, each observation contributes to the calculation of the increase, where
$g\simeq\frac{1}{N}\tilde{K}\sum_{l=1}^{N}(M_{l})^{T}h$, $g^{x}$ contribution
from initial condition increment, $g^{f}$ contribution from surface forcing
increment, $g^{b}$ contribution from open boundary increment and $M_{l}\equiv
M(t_{0},t_{0}+l\Delta t)$.
Moreover, we consider the sensitivity of $I$ to variations $\delta y_{0}$ in
the observations, in particular, we consider $\delta y_{0}=d$.

Figure 10 shows the total number of observations from each observing platform
that were assimilated into the model during cycle, increments and contribution
of observations from each platform to

    * –

NL: nonlinear model;

    * –

TL: total impact;

    * –

IC: initial condition;

    * –

FC: surface forcing conditions;

    * –

BC: open boundary conditions.

Figure 13 shows curve represents the evolution of the ocean state vector $I$
in time based on the either the prior or posterior control vector. In the case
of the observation impact, the actual contribution of each observation,
$y^{o}_{i}$, to the change $\Delta I$ defined in (28) (blue curve vs. black
curve) due to data assimilation is revealed. Conversely, the observation
sensitivity quantifies the change $\Delta I$ that will occur in $I$ a (red
curve vs. blue curve) as a result of perturbations in the observations,
$\delta y=\textbf{d}$. In Figures 15, 16 and 17 we show the results obtained
by running related tests on SCoPE.

Figure 7: IS4D-Var Cost Function: case1, case2, case3 and case 4,
respectively, using roms_wc13_2hours.in input file.

Figure 8: IS4D-Var Cost Function: case1, case2, case3 and case 4,
respectively, using roms_wc13_daily.in input file.

Figure 9: RBL4D-Var Cost functions using roms_wc13 3_2hours.in (left) and
roms_wc13_daily.in (right) file input, $Nouter=1$ and $Ninner=25$.

Figure 10: RBL4D-Var Cost functions using roms_wc13 3_2hours.in (left) and
roms_wc13_daily.in (right) file input, $Nouter=1$ and $Ninner=50$.

Figure 11: RBL4D-Var observation impact using roms_wc13_2hours.in (left) and
roms_wc13_daily.in (right) file. Figure 12: The 37N section along which the
time averaged transport is computed from the surface to a depth of 500m.
Figure 13: A schematic illustrating difference between the observation impact
and observation sensitivity calculations. Fixed time interval
$[t_{0},t_{0+7}]$, $y^{o}_{i}$ are observations at various times that it are
shown as pluses (+) while perturbed observations, $y^{o}_{i}+\delta y^{o}_{i}$
are indicated by filled circles ($\bullet$). Consequently, red curve is
perturbed analysis and blue curve is 4DVar analysis.

Figure 14: RBL4D-Var observation sensitivity using roms_wc13_2hours.in (left)
and roms_wc13_daily.in (right) file.

Figure 15: RBL4D-Var forecast impact using roms_wc13_2hours.in (left) and
roms_wc13_daily.in (right) file.

Figure 16: RBL4D-Var forecast impact in observations space using
roms_wc13_2hours.in (left) and roms_wc13_daily.in (right) file.

Figure 17: RBL4D-Var forecast sensitivity using roms_wc13_2hours.in (left) and
roms_wc13_daily.in (right) file.

##  4 Future developments

Main modification to apply DD in time in ROMS is described by point (c) in
Action 1.1, i.e. introduction of MPI communications in time. More precisely,
we need to introduce and manage the MPI communication in space and time. The
introduction of DD in time involves the following communications among
processes:

  * •

Intra communications: by splitting MPI communicator (OCN_COMM_WORLD) to create
local communicators (TASK_COMM_WORLD) to allow communications among processes
related to same spatial subdomain but different time intervals.
MPI commands are

    * –

MPI_Comm_split: partitions the group of MPI processes into disjoint subgroups
and creates a new communicator (TASK_COMM_WORLD) for each subgroup.

    * –

MPI_Isend and MPI_Irecv: sends and receives initial conditions by setting the
new communicator obtained from MPI_Comm_split.

  * •

Inter communications: by creating new communicators (OCNi_COMM_WORLD) to allow
communications among processes related to different spatial subdomains but
same time interval.
MPI commands are

    * –

MPI_Intercomm_create: creates an intercommunicator for each subgroup.

    * –

MPI_Isend and MPI_Irecv: sends and receives boundary conditions by setting the
new communicator obtained from MPI_Intercomm_create.

## References

  * [1] D’Amore, L., Cacciapuoti, R., 2021: “Model Reduction in Space and Time for decomposing ab initio 4D Variational Data Assimilation Problems”, Applied Numerical Mathematics, 2021, Volume 160, pp. 242-264 Elsevier, https://doi.org/10.1016/j.apnum.2020.10.003.
  * [2] https://www.myroms.org
  * [3] Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, and J. Zavala-Garay, 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems, Part I - System overview and formulation, Progress in Oceanography, 91, 34-49.
  * [4] Hedström, K.S., 2018: Technical Manual for a Coupled Sea-Ice/Ocean Circulation Model (Version 5), OCS Study BOEM 2018-007. U.S. Department of the Interior Bureau of Ocean Energy Management Alaska OCS Region. 182 p.
  * [5] Arango, H.G., Moore, A.M., Miller, A.J., Cornuelle, B.D., Lorenzo, E.Di, Cornuelle B.D. and Neilson, D.J., 2003: The ROMS Tangent Linear and Adjoint Models: A comprehensive ocean prediction and analysis system.
  * [6] D’Amore, L., Arcucci, R., Pistoia, J., Toumi, R. Murli, A. 2017: On the variational data assimilation problem solving and sensitivity analysis, Journal of Computational Physics, pp. 311-326.
  * [7] Langland, R. H., and N. Baker, 2004: Estimation of observation impact using the NRL atmospheric variational data assimilation adjoint system. Tellus, 56A, 189–201.
  * [8] Gürol, S., Weaver, A. T., Moore, A. M., Piacentini, A., Arango, and H. G., Gratton, S. 2014: B‐preconditioned minimization algorithms for variational data assimilation with the dual formulation. Quarterly Journal of the Royal Meteorological Society, 140(679), 539-556.


# Surface-directed and bulk spinodal decomposition compete to decide the
morphology of bimetallic nanoparticles

P. Pankaj  Saswata Bhattacharyya  Subhradeep Chatterjee
<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Materials
Science and Metallurgical Engineering,
Indian Institute of Technology, Hyderabad, 502284, Telangana, India

###### Abstract

An embedded-domain phase-field formalism is used for studying phase
transformation pathways in bimetallic nanoparticles (BNPs). Competition of
bulk and surface-directed spinodal decomposition processes and their interplay
with capillarity are identified as the main determinants of BNP morphology.
The former is characterized by an effective bulk driving force
$\Delta\tilde{f}$ which increases with decreasing temperature, while the
latter manifests itself through a balance of interfacial energies captured by
the contact angle $\theta$. The simulated morphologies, namely, core-shell,
Janus and inverse core-shell, cluster into distinct regions of the
$\Delta\tilde{f}$-$\theta$ space. Variation of $\theta$ with $\Delta\tilde{f}$
in the Ag-Cu alloy system is computed as a function of temperature using a
CALPHAD approach in which surface energies are estimated from a modified
Butler equation. This $\theta-\Delta\tilde{f}$ trajectory for Ag-Cu, when
superimposed on the morphology map, enables the prediction of different
morphological transitions as a function of temperature. Therefore, the study
establishes a unique thermodynamic framework coupled with phase-field
simulations for predicting and tailoring nanoparticle morphology through a
variation of processing temperature.

## Introduction

Properties of bimetallic nanoparticles (BNPs) used in diverse fields such as
catalysis, photonics, spintronics and biomedical-sensing [1, 2, 3] depend
crucially on their morphology. Core-shell (CS) and Janus are the most commonly
reported morphologies: CS morphology consists of an outer shell of the lower
surface energy component surrounding a core of the other component, while
Janus is characterized by the two components forming two opposite faces of a
particle with their common interface extending to the surface. An inverse
core-shell (ICS) morphology, where the higher surface energy component forms
the shell, has also been observed in a few cases [4, 5].

A fundamental understanding of morphological development in BNPs is crucial
for tailoring their properties. First principles [6] and atomistic [7]
simulations, as well as classical thermodynamics [8] have been employed to
find the lowest energy morphology of BNPs. However, these approaches generally
do not address kinetic aspects of morphological development. Since many
commonly used BNP systems like Ag-Cu, Ag-Ni, Au-Co, Co-Cu, etc. exhibit solid-
state immiscibility, spinodal decomposition (SD) presents a likely kinetic
pathway for morphology development in BNPs [9]. Recently, an embedded-domain
phase-field model (EPFM) was employed [10] to study morphological development
in BNPs through SD and understand how contact angle and particle size
influence it.

From a processing perspective, temperature, or equivalently, undercooling
below a critical temperature, is a variable that exerts a strong influence on
the phase separation process by altering the driving force for bulk SD. In
addition, in confined systems such as BNPs, surface directed SD (SDSD)
presents another mode of microstructural evolution. Interplay between these
two alternative mechanisms and their interactions with capillarity give rise
to different BNP morphologies. In this paper, we use EPFM to systematically
investigate this process that offers further insight into morphology selection
in BNPs.

## Model Formulation

In the EPFM formalism [10] illustrated schematically in Fig. 1, a BNP system
consisting of atomic species A and B is modeled by placing an isolated
particle ($\beta$) in an inert matrix ($\alpha$). A scaled composition field
$c$ describes phase separation of an initially homogeneous $\beta$ to a
mixture of $\beta_{1}$ and $\beta_{2}$ phases with equilibrium compositions
$c_{\beta_{1}}^{e}$ and $c_{\beta_{2}}^{e}$, respectively. We use
$c=(X-X_{\beta_{1}}^{e})/(X_{\beta_{2}}^{e}-X_{\beta_{1}}^{e})$ for the
scaling where $X$ denotes the mole fraction of B and $X_{\beta_{i}}$ is the
equilibrium composition of phase $\beta_{i}$ at a given temperature. An
auxiliary non-conserved phase-field variable $\phi$ is used to distinguish the
particle from the matrix; we use a stationary and radially symmetric tanh
profile for $\phi$ that varies smoothly across the particle surface from 0 in
the matrix side to 1 inside the particle [10]. Total free energy of the system
is then expressed as:

$\mathcal{F}=\frac{1}{V_{m}}\int_{\Gamma}\Big{(}f(c,\phi)+\kappa_{\phi}|\nabla\phi|^{2}+\kappa_{c}|\nabla
c|^{2}\Big{)}d\Gamma,$ (1)

where $V_{m}$ is the molar volume, $f$ is bulk free energy density,
$\kappa_{c}$ and $\kappa_{\phi}$ are gradient energy coefficients associated
with $c$ and $\phi$, respectively, and $\Gamma$ is the volume of the whole
computational domain.

Figure 1: Computational domain $\Gamma$ consisting of an isolated nanoparticle
$\beta$ (physical domain of interest) in an inert matrix $\alpha$. Periodic
boundary conditions are imposed on the domain boundary $d\Gamma$.

Defining an interpolating function $h(\phi)=\phi^{3}(10-15\phi+6\phi^{2})$
that connects $\alpha$ and $\beta$ (with $\phi$ being 0 and 1 inside $\alpha$
and $\beta$, respectively), the bulk free energy density $f$ is given as
$f(c,\phi)=h(\phi)f^{\beta}+(1-h(\phi))f^{\alpha}+\omega(c)g(\phi)$. Here
$\omega(c)g(\phi)=\omega_{0}(1-\chi c)\phi^{2}(1-\phi)^{2}$ describes the free
energy barrier between $\alpha$ and $\beta$, $\omega_{0}$ sets the barrier
height and $\chi$ controls solute segregation at the particle-matrix interface
[10]. We use the following forms of free energy for $\alpha$ and $\beta$:

$\displaystyle f^{\alpha}$ $\displaystyle=f_{0}^{m}(c-c_{\alpha}^{e})^{2},$
$\displaystyle f^{\beta}$
$\displaystyle=f_{0}^{p}(c-c_{\beta_{1}}^{e})^{2}(c-c_{\beta_{2}}^{e})^{2}.$
(2)

The scaled equilibrium compositions
$c_{\alpha}^{e},c_{\beta_{1}}^{e},c_{\beta_{2}}^{e}$ are chosen to be 0.5, 0
and 1, respectively, and $f_{0}^{p}$ and $f_{0}^{m}$ are temperature-dependent
constants. Evolution of composition field within the particle is described by
the Cahn-Hilliard equation:

$\frac{\partial{c}}{\partial{t}}=\nabla\cdot{M(\phi)}\nabla\frac{\delta\mathcal{F}}{\delta{c}}.$
(3)

We constrain the matrix to remain inert with respect to solute diffusion by
choosing $M(\phi)=M_{c}h(\phi)$, with $M_{c}$ being the effective atomic
mobility of solute in the particle. Eq. (3) is solved numerically using a
semi-implicit Fourier spectral method [11, 12, 13]. We non-dimensionalize all
parameters in the equations using characteristic length $L_{c}$, time
$\tau_{c}$ and energy $E_{c}$. Details of the non-dimensionalization procedure
are provided Supplementary Information A; further details of the model and its
numerical implementation can be found elsewhere [10].

Energies of $\alpha$-$\beta_{1}$, $\alpha$-$\beta_{2}$ and
$\beta_{1}$-$\beta_{2}$ interfaces, designated as $\sigma_{1},\sigma_{2}$ and
$\sigma_{12}$, are obtained from the equilibrium composition profiles
$c^{e}(x)$ across the respective interfaces:

$\displaystyle\sigma$
$\displaystyle=\frac{1}{V_{m}}\int\Big{[}{f}(c^{e}(x),\phi(x))+\kappa_{\phi}|\nabla\phi|^{2}+\kappa_{c}|\nabla{c^{e}(x)}|^{2}$
$\displaystyle-\Big{\\{}(1-c^{e}(x))\mu_{A}^{e}+c^{e}(x)\mu_{B}^{e}\Big{\\}}\Big{]}dx,$
(4)

where $\mu_{i}^{e}\;(i=A,\;B)$ denotes the equilibrium chemical potential of
component $i$ in any one the coexisting phases across the interface.

The surface and interfacial energies define the contact angle $\theta$ at the
triple junction between the phases as:

$\cos\theta=\frac{\sigma_{1}-\sigma_{2}}{\sigma_{12}}.$ (5)

Phase transformations within the particle need not always yield triple
junctions at the surface. Cahn, in his classic paper on “Critical Point
Wetting”, described a condition that precluded the formation of a triple
junction. He defined this _spontaneous_ wetting condition to be
$\sigma_{1}-\sigma_{2}\geq\sigma_{12}$; the equality condition corresponds to
$\theta=$$$, while $\theta$ become undefined for the inequality condition. In
the latter case, the phase with the low surface energy develops a continuous
layer on the surface.

Non-zero $\chi$ leads to preferential solute segregation to one of the
surfaces by creating imbalance between the ‘surface’ energies $\sigma_{1}$ and
$\sigma_{2}$ – its key role on BNP morphology has already been explored in
detail [10]. Here, we keep $\chi$ and particle size fixed, and investigate the
competition of bulk SD and SDSD in deciding the BNP morphology (CS, ICS and
Janus) by varying $f_{0}^{p}$, $\omega_{0}$ and $\kappa_{c}$. Note that these
variable also control $\theta$. Since each of these parameter combinations
gives rise to one of these three morphologies, we performed simulations over a
large set of parameters and classified the results in terms of driving force
and contact angle.

## Results and Discussion

###  .1 Simulation of BNP morphologies

All simulations begin with an initially homogeneous, axisymmetric $\beta$
particle with composition $c=0.5$ quenched inside the miscibility gap. A small
initial noise ($\pm{1}\%$) mimicking thermal fluctuations is applied to
trigger bulk SD. Figs. 2 and 3 present time snapshots of evolution of stable
and metastable morphologies corresponding to representative parameter sets
listed in Table 1. In these figures, green color represents the solute-poor
$\beta_{1}$ phase, blue the solute-rich $\beta_{2}$ phase, and grey the
undecomposed $\beta$ with $c=0.5$.

Table 1: Variation of surface and interfacial energies and contact angle $\theta$ with model parameters. Figure numbers refer to snapshots of morphological evolution for the corresponding set of parameters. Set | $f_{0}^{p}$ | $\kappa_{c}$ | $\omega_{0}$ | $\sigma_{1}$ | $\sigma_{2}$ | $\sigma_{12}$ | $\theta$ | Figure no.
---|---|---|---|---|---|---|---|---
1 | 8 | 1 | 12 | 4.19 | 3.09 | 0.94 | wetting |  2(a)-2(e)
2 | 6 | 2 | 3.75 | 1.77 | 1.44 | 1.15 | $None$ |  2(f)-2(j)
3 | 4 | 1 | 6 | 2.21 | 1.76 | 0.67 | $None$ |  3(a)-3(e)
4 | 2 | 8 | 6 | 2.24 | 1.86 | 1.33 | $None$ |  3(f)-3(j)
5 | 4 | 2 | 5 | 1.98 | 1.60 | 0.94 | $None$ |  3(k)-3(o)

In all cases, SDSD precedes bulk-SD and forms alternate solute-rich and
solute-lean rings which grow inward. When bulk SD starts in the interior of
the particle, it creates intertwined compositionally modulated domains that
interact with the rings growing from the surface to the center. Depending on
the chosen parameter set, domain coarsening proceeds along different pathways
to produce the final BNP configuration. A Janus structure is formed when
coarsening disrupts the continuity of ring-like structures at the surface, and
CS or ICS results otherwise.

Figs. 2(a-e) is an example of a typical evolution pattern leading to a stable
CS morphology (set 1 of Table 1). It occurs irrespective of the driving force
if the spontaneous wetting condition is satisfied (_i.e._ ,
$\sigma_{1}-\sigma_{2}\geq\sigma_{12}$). For large $\theta$, if bulk driving
force is sufficiently high to break the outermost layer, we obtain stable
Janus morphology. A typical sequence of its formation is presented in Figs.
2(f-j) which corresponds to parameter set 2 of Table 1.

(a) $t=0.2$

(b) $t=1$

(c) $t=20$

(d) $t=600$

(e) $t=3500$

(f) $t=0.4$

(g) $t=4$

(h) $t=200$

(i) $t=2000$

(j) $t=$3.9E4

Figure 2: Time snapshots of morphological evolution of stable CS (top row) and
Janus (bottom row) configurations. Non-dimensional times are indicated below
each snapshot.

Parameter sets 3-5 of Table 1 correspond to intermediate values of driving
force (which scales with $f_{0}^{p}$) and $\theta$ which result in metastable
configurations shown in Figure 3. For example, Figs. 3(a-e) shows how a CS
morphology can develop for small non-zero $\theta$’s if the driving force for
bulk spinodal is also small. Thus, the spontaneous wetting condition is only a
sufficient condition – stable CS forms when it is satisfied, but a metastable
CS can develop even when it is not.

(a) $t=0.4$

(b) $t=4$

(c) $t=200$

(d) $t=2000$

(e) $t=6000$

(f) $t=2$

(g) $t=12$

(h) $t=80$

(i) $t=600$

(j) $t=19000$

(k) $t=1$

(l) $t=10$

(m) $t=200$

(n) $t=2000$

(o) $t=55000$

Figure 3: Time snapshots of metastable CS (top row) and ICS (middle and bottom
rows) evolution. Non-dimensional times are indicated below for each snapshot.

On the other hand, one obtains a metastable ICS morphology when the driving
force is small but $\theta$ is large, as demonstrated by the typical evolution
patterns in Figs. 3(f-j) and Figs. 3(k-o). When surface effects dominate over
the bulk, an onion-like ring structure develops initially (Figs. 3(a)-3(b)),
which coarsens subsequently (Figs. 3(a-e)) to form a metastable ICS
morphology. It is also possible that we get metastable ICS/CS even in the
presence of bulk spinodal. As shown in Figs. 3(f-o), the bulk driving force in
this case is not sufficiently high enough to achieve a stable Janus
configuration, and domain coarsening leads to an ICS configuration. Between
the two evolutionary paths to ICS, the latter is observed when the driving
force is higher.

As shown in Fig. 4, however, the CS at non-zero $\theta$ and ICS
configurations are metastable, as they relax to more stable structures when
subjected to a sustained white noise [14]. Here the initial configurations for
the simulations are the final metastable ones shown in Fig. 3. With increasing
time, the concentric structures break down and evolve to form Janus. The time
taken for the metastable-to-stable morphological transition is very large
because the difference of the energies between the metastable (CS/ICS) and
stable (Janus) configuration is very small, resulting in very sluggish
diffusion. Therefore, these can be termed as kinetically trapped
configurations [15, 16, 10].

(a) $t=6000$

(b) $t=6300$

(c) $t=7000$

(d) $t=10000$

(e) $t=20000$

(f) $t=19000$

(g) $t=20600$

(h) $t=25000$

(i) $t=40000$

(j) $t=111600$

(k) $t=55000$

(l) $t=57000$

(m) $t=70000$

(n) $t=82000$

(o) $t=127000$

Figure 4: Time snapshots of transition from metastable CS/ICS to stable CS
configurations brought about by sustained noise. Non-dimensional times are
indicated below for each snapshot.

To understand the formation of different morphologies as a function of the
model parameter sets, one needs to relate the competing processes of bulk SD
and SDSD. For bulk SD, $f_{0}^{p}$ sets the driving force $\Delta{f}$, defined
as the difference between the free energy of non-equilibrium initial $\beta$
state and that of the final state made up of equal parts of equilibrium
$\beta_{1}$ and $\beta_{2}$ phases:

$\displaystyle\Delta{f}$
$\displaystyle=f^{\beta}(c=0.5)-\frac{1}{2}\left[f^{\beta}(c=c_{\beta_{1}}^{e})+f^{\beta}(c=c_{\beta_{2}}^{e})\right]$
$\displaystyle=f_{0}^{p}/16\,\,(\textrm{using Eq.~{}\eqref{eq:bulk-free-
energy-dens}}).$ (6)

$\Delta{f}$, together with $\kappa_{c}$, also determines the interfacial
energy $\sigma_{12}$ of the system as [17]:

$\displaystyle\sigma_{12}=\frac{1}{3V_{m}}\left[c_{\beta_{2}}^{e}(T)-c_{\beta_{1}}^{e}(T)\right]^{3}\sqrt{\kappa_{c}\Delta{f(T)}}.$
(7)

Both $\Delta{f}$ and $\sigma_{12}$ increase with decreasing $T$ (or increasing
undercooling $\Delta{T}$).

The surface energies $\sigma_{1}$ and $\sigma_{2}$, and their difference
$\Delta\sigma$, on the other hand, scale with the height of the free energy
barrier $\omega_{0}$. While $\Delta f$ controls bulk SD within the particle,
$\omega_{0}$ influences SDSD by defining $\Delta\sigma$ – the larger the
$\Delta\sigma$, the easier it is to initiate SDSD. We now define a normalized
driving force for bulk phase separation, $\Delta\tilde{f}$, by taking the
ratio of $\Delta f$ to $\Delta\sigma$.

Fig. 5 presents the essence of all simulations by assigning the location of
BNP configuration in the space of model parameters $\Delta\tilde{f}$ and
$\theta$. The map clearly demarcates three distinct regions in this space:
high $\Delta\tilde{f}$-high $\theta$ for Janus, low $\Delta\tilde{f}$-high
$\theta$ for ICS, and the intervening region for CS. The region of high
$\Delta\tilde{f}$-low $\theta$ is physically inaccessible since high
$\Delta\tilde{f}$ implies either high $\Delta{f}$ or low $\omega_{0}$, which
increases $\sigma_{12}$ or decreases $\Delta\sigma$, respectively, and thereby
lead to high $\theta$.

(a) Figure 5: Simulated stable and metastable BNP configurations in the
$\Delta\tilde{f}$–$\theta$ space. Dashed lines are drawn to delineate regions
of where Janus, CS and ICS morphologies are observed.

When the spontaneous wetting condition is satisfied ($\theta=$$$ or
undefined), CS is the final configuration irrespective of $\Delta\tilde{f}$.
Also, for a given $\theta$, a metastable CS configuration is preferred over
ICS at larger $\Delta\tilde{f}$. At lower $\Delta\tilde{f}$, SDSD creates
alternate concentric rings extending towards the center, with segregation
giving rise to an outermost $\beta_{1}$ ring [10]. Subsequently, ring-
coarsening leads to the ICS structure. In contrast, bulk SD dominates at
higher $\Delta\tilde{f}$, resulting in interactions between inner
interconnected domains and the outermost $\beta_{1}$ ring during coarsening.
Subsequently, the latter is replaced by $\beta_{2}$, leading to metastable CS.
Note that these interactions do not yield the stable Janus configuration,
because, in this case, $\Delta\sigma$ remains sufficiently large to cause
$\beta_{2}$ to preferentially spread along the surface and $\beta_{1}$ to
recede.

###  .2 Computation of chemical and capillary forces for Ag-Cu

Fig. 5 demonstrated how different morphologies clustered around in distinct
regions of the driving force – contact angle space. This space essentially
represents the interplay of chemical and capillary forces, and temperature is
one of the key physical parameters that directly or indirectly controls them.
The chemical driving force for phase separation, $\Delta{f}$, for a given
system can be obtained as a function of temperature in a relatively
straightforward way from its CALPAHD data. Capillarity, on the other hand, is
manifested through the contact angle $\theta$, which itself is defined in
terms of interfacial and surface energies
($\sigma_{1},\sigma_{2},\sigma_{12}$). Reliable values of the latter as a
function of temperature are often difficult to measure experimentally.
Nevertheless, here we attempt to arrive at fair estimates of these energies
using the available data and thermodynamic correlations. We consider the Ag-Cu
system as an example and proceed to compute the temperature dependence of
these quantities.

####  .2.1 Driving force for phase separation

The molar bulk free energy $F_{m}$ of the Cu-Ag solid solution can be
expressed through a Redlich-Kister polynomial as:

$\displaystyle F_{m}=$
$\displaystyle\left[X_{\mathrm{Ag}}F_{\mathrm{Ag}}^{0}+X_{\mathrm{Cu}}F_{\mathrm{Cu}}^{0}\right]+RT\left[X_{\mathrm{Cu}}\ln{X_{\mathrm{Cu}}}+X_{\mathrm{Ag}}\ln{X_{\mathrm{Ag}}}\right]$
$\displaystyle+X_{\mathrm{Cu}}X_{\mathrm{Ag}}\left[L_{0}+L_{1}(X_{\mathrm{Ag}}-X_{\mathrm{Cu}})\right]$
(8)

where $X_{i}$’s ($i=$ Cu, Ag) are the mole fractions component $i$ in the
solution. The first bracketed group in the right hand side is the mixture free
energy contribution with $F_{i}^{0}$’s being the standard free energies of the
pure components which can be found in the SGTE database [18]. The second group
is the contribution from the ideal solution part of free energy
($F_{m}^{\mathrm{id,bulk}}$) and the last term represents the excess
contribution ($F_{m}^{\mathrm{ex,bulk}}$). The temperature-dependent
interaction parameters for the Cu-Ag solid solution are given as
$L_{0}=34532-9.178T$ and $L_{1}=-5996+1.725T$ [19].

Equilibrium solute mole fractions $X_{i}^{e}$ in Cu-rich and Ag-rich solid
solutions below the critical temperature $T_{c}$ are computed using the
conditions for chemical equilibrium:

$\displaystyle\mu_{\mathrm{Ag}}^{\mathrm{Cu_{ss}}}=\mu_{\mathrm{Ag}}^{\mathrm{Ag_{ss}}},\quad\mu_{\mathrm{Cu}}^{\mathrm{Cu_{ss}}}=\mu_{\mathrm{Cu}}^{\mathrm{Ag_{ss}}},$
(9)

where the chemical potentials $\mu_{i}\,(i=\mathrm{Cu},\mathrm{Ag})$ evaluated
at equilibrium bulk compositions $X_{i}^{e}$ are given by

$\displaystyle\mu_{\mathrm{Cu}}=F_{m}-X_{\mathrm{Ag}}\frac{\partial
F_{m}}{\partial
X_{\mathrm{Ag}}};\quad\mu_{Ag}=F_{m}+(1-X_{\mathrm{Ag}})\frac{\partial
F_{m}}{\partial X_{\mathrm{Ag}}}.$ (10)

The free energy of mixing for an alloy composition $X_{\mathrm{Ag}}$ is
expressed as

$\displaystyle\Delta{F}_{\mathrm{mix}}=F_{m}-(X_{\mathrm{Cu}}\mu_{\mathrm{Cu}}^{e}+X_{\mathrm{Ag}}\mu_{\mathrm{Ag}}^{e}).$
(11)

The driving force for bulk spinodal, $\Delta{f}$, is the maximum value of
$\Delta{F}_{\mathrm{mix}}$ over the entire composition range (_i.e._ , maximum
of the $\Delta{F}_{\mathrm{mix}}-X$ curve). It is evaluated and plotted in
Fig. 6(a) for the temperature range of 400-800 K. The plot shows that
$\Delta{f}$ decreases monotonically with increasing temperature. Using the
same thermodynamic data, we also compute the equilibrium mole fractions
$X_{i}^{e}$ ($i=\;\mathrm{Cu,Ag}$) in both Cu- and Ag-rich solutions which are
required for estimating the interfacial energy.

####  .2.2 Interfacial energy

The temperature dependence of the interfacial energy between the Cu- and Ag-
rich solid solutions ($\sigma_{\mathrm{Cu-Ag}}$) can be computed using Eq.
(7). To be able to do so, however, one requires the value of $\kappa_{c}$. We
utilize the value of $\sigma_{\mathrm{Cu-Ag}}=0.197$
$\mathrm{J}\mathrm{m}^{-2}$ at 800 K obtained from molecular dynamics
simulations [7], along with the values of $\Delta f$ and $X_{i}^{e}$ at this
temperature calculated in the previous step, and a constant molar volume of
$10^{-5}$ $\mathrm{m}^{3}\mathrm{mol}^{-1}$, to compute $\kappa_{c}=1.94\times
10^{-14}$ $\mathrm{J}\mathrm{m}^{2}\mathrm{mol}^{-1}$. Since $\kappa_{c}$ is
generally considered to be temperature-independent for a given system, we use
this constant value in Eq. (7) to compute $\sigma_{\mathrm{Cu-Ag}}$ in the
temperature range of 400-1200 K.

(a)

(b)

Figure 6: (a) Variation of $\Delta{f}$ with temperature for the Ag-Cu system
computed using its CALPHAD data. (b) Temperature-dependence of interfacial and
surface energies of Cu-poor and Cu-rich solid solution phases. Temperature of
transition ($T_{tr}$) from stable CS to Janus morphology is indicated.

####  .2.3 Surface energies

We calculate surface energies $\sigma_{i,\mathrm{ss}}$ of the Cu-rich and Ag-
rich solid solutions as a function of temperature $T$ using a modified Butler
model [20]:

$\displaystyle\sigma(T)=$
$\displaystyle\sigma_{i}^{0}+\frac{RT}{A_{i}}\log{\left[\frac{X_{i}^{\mathrm{surf}}}{X_{i}^{e}}\right]}$
$\displaystyle+$
$\displaystyle\frac{1}{A_{i}}[F_{i}^{\mathrm{ex,surf}}(T,X_{i}^{\mathrm{surf}})-F_{i}^{\mathrm{ex,bulk}}(T,X_{i}^{e})].$
(12)

Here, $i=\mathrm{A,\;B}$ denote the components Cu and Ag, respectively,
$\sigma_{i}^{0}$ are the surface energies of the pure components, $A_{i}$ are
the molar surface areas, $X_{i}^{\mathrm{surf}}$ are the surface compositions
at $T$, $F_{i}^{\mathrm{ex,surf}}$ and $F_{i}^{\mathrm{ex,bulk}}$ are the
excess partial molar free energies of $i$ associated with surface and bulk,
respectively.

Surface energies of the pure components Eq. (12) are obtained from the
correlation

$\displaystyle\sigma_{i}^{0}=1.25\sigma_{i}^{\mathrm{liq}}+\frac{d\sigma}{dT}(T-T_{m}^{i})$
(13)

where $\sigma_{i}^{\mathrm{liq}}$ is the surface tension of liquid component
at its melting point $T_{m}^{i}$ and the temperature coefficient of surface
energy is taken as $-10^{-4}$ $\mathrm{J}\mathrm{m}^{-2}\mathrm{K}$. The
partial molar excess free energies for bulk, $F_{i}^{\mathrm{ex,bulk}}$, can
be expressed in terms of $X_{i}^{e}$ using Eq. (10), with excess free energy
$F^{ex,bulk}$ replacing the total molar free energy $F_{m}$. Following Tanaka
and Hara [20], we take
$F_{i}^{\mathrm{ex,surf}}=\beta_{mix}F_{i}^{\mathrm{ex,bulk}}$. The parameter
$\beta_{mix}$ is the ratio of coordination number in the surface to that in
the bulk; it is taken as $0.75$ for face centered cubic solid solutions. Molar
surface areas $A_{i}$ of the pure components are obtained using the relation
$A_{i}=1.091N_{0}^{1/3}V_{m}^{2/3}$ where $N_{0}$ is the Avogadro’s number.
Eq. (12) constitutes a set of two simultaneous non-linear algebraic equations
which can now be solved numerically to obtain the unknown surface composition
$X_{i}^{\mathrm{surf}}$ for both the terminal phases. Plugging its values in
Eq. (10), one can compute the surface energies of the Cu- and Ag-rich solution
phases.

Fig. 6(b) plots the difference of the surface energies,
$\Delta\sigma=\sigma_{\mathrm{Cu_{ss}}}-\sigma_{\mathrm{Ag_{ss}}}$, and the
interfacial energy, $\sigma_{\mathrm{Cu-Ag}}$, with temperature. It shows that
although the surface energies themselves vary with temperature, their
difference is relatively insensitive to it, changing only slightly at higher
temperatures. The interfacial energy, on the other hand, depends strongly with
temperature, decreasing steeply at higher temperatures. This is expected, as
it must vanish at the critical temperature for the miscibility gap. In terms
of the spontaneous wetting criterion (Eq. (5)), temperature of intersection of
surface and interfacial energy lines, $T_{tr}$ ($\sim$ 910 K) marks the
transition from _stable_ CS to Janus morphologies.

###  .3 Connecting phase-field results with thermodynamic computations

We proceed further to obtain a correlation between the effective driving force
$\Delta\tilde{f}$ and $\theta$, noting that the latter is defined only for
$T\leq T_{tr}$. First, the bulk driving force $\Delta{f}$ (in J/mol) for phase
separation in Ag-Cu is normalized by the product of surface energy difference
$\Delta\sigma$ (in J/m2) and molar surface area (in m2/mol). Fig. 7(a) shows
the temperature dependence of $\Delta\tilde{f}$ estimated this way, it is very
similar to the variation of $\Delta{f}$ with $T$ shown in Fig. 6(a). Next, the
interfacial and surface energies computed earlier are used to determine how
contact angle varies as a function of temperature. This is presented in Fig.
7(b), which shows that at high temperatures, $\theta$ decreases steeply with
decrease in temperature, but the rate of this decrease reduces at lower
temperatures. Finally, these variations are combined into a single
$\Delta\tilde{f}$-$\theta$ plot in Fig. 7(c). It shows a monotonic increase of
$\theta$ with $\Delta\tilde{f}$; however, the slope of the curve is steeper at
low-$\Delta\tilde{f}$, gradually becoming gentler with increase in
$\Delta\tilde{f}$ (_i.e._ , decrease in $T$). This line is the trajectory that
the system follows in the $\Delta\tilde{f}$-$\theta$ space as temperature is
reduced, and therefore captures its response to a change in the state variable
$T$. This trajectory is now superimposed on the morphology map which is
redrawn in Fig. 7(d). It predicts the morphological transitions in Ag-Cu for
$T<T_{tr}$: both $\Delta\tilde{f}$ and $\theta$ increase with decreasing
temperature, and the system moves from a metastable CS to metastable ICS
configuration, before finally transitioning into the Janus regime. Thus,
depending on the processing conditions, all the three configurations can form
in Ag-Cu BNPs, as confirmed by the experimentally observed configurations [5,
21].

Since the temperature dependence of bulk chemical and capillary forces are
system specific, it should be noted that the trajectory may not always pass
through the metastable ICS region of the map for all alloy systems. However,
the insights gained from the study remains valid and provide crucial
guidelines and understanding for conducting further experiments aimed at
tailoring the BNP morphology.

(a)

(b)

(c)

(d)

Figure 7: (a) Temperature-dependence of $\Delta{\tilde{f}}$. (b) Variation of
$\theta$ with temperature for the Ag-Cu system computed using its CALPHAD
data. (c) Correlation of $\theta$ with $\Delta{\tilde{f}}$. (d) The
$\theta-\Delta\tilde{f}$ trajectory for Ag-Cu is superimposed on the
morphology map obtained from simulations. CS$\to$ICS and ICS$\to$Janus
transitions are indicated by the crosses, and vertical dashed lines are drawn
on (a) and (b) to indicate the temperatures for these transitions (the former
takes place at the higher $T$).

##  I Conclusions

Results of the present study show how different BNP morphologies can emerge
from a competition between two alternative and concomitant mechanisms, namely,
bulk and surface-directed spinodal decomposition. Their interplay with
capillarity sets up the eventual coarsening pathway toward a steady-state
configuration. When we express the results in terms of three physical
parameters, namely, driving force, difference in surface energies and contact
angle, different morphologies automatically cluster into three distinct
regions in the $\Delta\tilde{f}-\theta$ space. This identification of relevant
physical parameters appears remarkable, as the computed thermodynamic
trajectory for Ag-Cu involving the same variables traverses all the three
distinct regions revealing the morphological transitions. We note that
morphological transitions for a particular alloy system are sensitive to the
nature of temperature-dependence of the relevant variables. Therefore, the
exact transition points for different alloys systems will be different.

The following specific conclusions can be drawn from the study:

  1. 1.

Irrespective of the driving force for bulk spinodal, stable CS forms when the
spontaneous wetting condition is satisfied ($\theta=$$$). This happens at
$T\geq{T_{tr}}$.

  2. 2.

A combination of low driving force and high contact angle gives rise to
metastable ICS, while metastable CS forms at moderate driving force and lower
non-zero $\theta$. Thus, the _spontaneous_ wetting condition given by Cahn is
found to be a _sufficient_ , but not a _necessary_ condition for the formation
of CS.

  3. 3.

Janus forms when both bulk driving force and contact angle are large.

  4. 4.

Trajectory of Ag-Cu in the driving force-contact angle space shows transitions
involving CS, ICS and Janus morphologies as a function of temperature.

## Conflict of Interest

The authors have no conflicts to disclose.

## Acknowledgement

Authors gratefully acknowledge the computational support from DST-NSM Grant
DST/NSM/R&D-HPC-Applications/2021/03.

## References

  * Ferrando [2016] R. Ferrando,  _Structure and Properties of Nanoalloys_ (Elsevier, 2016).
  * Cortie and McDonagh [2011] M. B. Cortie and A. M. McDonagh, Synthesis and optical properties of hybrid and alloy plasmonic nanoparticles, Chemical Reviews 111, 3713 (2011).
  * Medina-Cruz  _et al._ [2020] D. Medina-Cruz, B. Saleh, A. Vernet-Crua, A. Nieto-Argüello, D. Lomelí-Marroquín, L. Y. Vélez-Escamilla, J. L. Cholula-Díaz, J. M. García-Martín, and T. Webster, Bimetallic nanoparticles for biomedical applications: A review, in  _Racing for the Surface: Antimicrobial and Interface Tissue Engineering_ , edited by B. Li, T. F. Moriarty, T. Webster, and M. Xing (Springer International Publishing, Cham, 2020) pp. 397–434.
  * Tsuji  _et al._ [2010] M. Tsuji, S. Hikino, R. Tanabe, and D. Yamaguchi, Synthesis of Ag@Cu core–shell nanoparticles in high yield using a polyol method, Chemistry Letters 39, 334 (2010).
  * Nakamura  _et al._ [2007] T. Nakamura, Y. Tsukahara, T. Yamauchi, T. Sakata, H. Mori, and Y. Wada, Preparation of Ag core–Cu shell nanoparticles by microwave-assisted alcohol reduction process, Chemistry Letters 36, 154 (2007).
  * Ferrando [2015] R. Ferrando, Symmetry breaking and morphological instabilities in core-shell metallic nanoparticles, Journal of Physics Condensed Matter 27, 13003 (2015).
  * Chandross [2014] M. Chandross, Energetics of the formation of Cu–Ag core–shell nanoparticles, Modelling and Simulation in Materials Science and Engineering 22, 075012 (2014).
  * Yuan  _et al._ [2008] C. W. Yuan, S. J. Shin, C. Y. Liao, J. Guzman, P. R. Stone, M. Watanabe, J. W. Ager, E. E. Haller, and D. C. Chrzan, Structure map for embedded binary alloy nanocrystals, Applied Physics Letters 93, 193114 (2008).
  * Radnóczi  _et al._ [2017] G. Radnóczi, E. Bokányi, Z. Erdélyi, and F. Misják, Size dependent spinodal decomposition in Cu-Ag nanoparticles, Acta Materialia 123, 82 (2017).
  * Pankaj  _et al._ [2022] P. Pankaj, S. Bhattacharyya, and S. Chatterjee, Competition of core-shell and Janus morphology in bimetallic nanoparticles: Insights from a phase-field model, Acta Materialia 233, 117933 (2022).
  * Cogswell and Carter [2011] D. A. Cogswell and W. C. Carter, Thermodynamic phase-field model for microstructure with multiple components and phases: The possibility of metastable phases, Physical Review E 83, 061602 (2011).
  * Cogswell [2010] D. A. Cogswell,  _A phase-field study of ternary multiphase microstructures_ , Ph.D. thesis, Massachusetts Institute of Technology (2010).
  * Zhu  _et al._ [1999] J. Zhu, L.-Q. Chen, J. Shen, and V. Tikare, Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: Application of a semi-implicit Fourier spectral method, Physical Review E 60, 3564 (1999).
  * Cook [1970] H. Cook, Brownian motion in spinodal decomposition, Acta Metallurgica 18, 297 (1970).
  * Li  _et al._ [2019] L. Li, X. Li, Z. Duan, R. J. Meyer, R. Carr, S. Raman, L. Koziol, and G. Henkelman, Adaptive Kinetic Monte Carlo simulations of surface segregation in PdAu nanoparticles, Nanoscale 11, 10524 (2019).
  * Grammatikopoulos  _et al._ [2016] P. Grammatikopoulos, J. Kioseoglou, A. Galea, J. Vernieres, M. Benelmekki, R. E. Diaz, and M. Sowwan, Kinetic trapping through coalescence and the formation of patterned Ag–Cu nanoparticles, Nanoscale 8, 9780 (2016).
  * Haider  _et al._ [2007] F. Haider, R. Kozubski, and T. Abinandanan, Simulation techniques, in  _Alloy Physics_ (John Wiley & Sons, Ltd, 2007) Chap. 12, pp. 653–706.
  * Dinsdale [1991] A. Dinsdale, Sgte data for pure elements, Calphad 15, 317 (1991).
  * Subramanian and Perepezko [1993] P. Subramanian and J. Perepezko, The Ag-Cu (silver-copper) system, Journal of Phase Equilibria 14, 62 (1993).
  * Tanaka and Hara [2001] T. Tanaka and S. Hara, Thermodynamic evolution of nanoparticle binary alloy phase diagrams, Zeitschrift fuer Metallkunde/Materials Research and Advanced Techniques 92, 1236 (2001).
  * Malviya and Chattopadhyay [2014] K. D. Malviya and K. Chattopadhyay, Synthesis and mechanism of composition and size dependent morphology selection in nanoparticles of Ag–Cu alloys processed by laser ablation under liquid medium, The Journal of Physical Chemistry C 118, 13228 (2014).
  * C.P. Wang  _et al._ [2012] C.P. Wang, L.N. Yan, J.J. Han, and X.J. Liu, Diffusion mobilities in the fcc Ag–Cu and Ag–Pd alloys, Calphad 37, 57 (2012).
  * Cahn [1961] J. W. Cahn, On spinodal decomposition, Acta Metallurgica 9, 795 (1961).
  * Cahn [1977] J. W. Cahn, Critical point wetting, The Journal of Chemical Physics 66, 3667 (1977).

## Appendix: Non-dimensionalization procedure

We make the governing equations (Eqs. (1) and (3)) dimensionless by using
dividing the length, energy and time variables by their characteristic values
($L_{c}$, $E_{c}$, and $\tau_{c}$, respectively):

$\displaystyle x=x^{\prime}/L_{c},f=f^{\prime}/E_{c},t=t^{\prime}/\tau_{c},$
(14)

where the primed quantities represent the dimensional values of the variables.
Now we write the dimensional form of Eq. (1) and use the above relations to
make it non-dimensional (noting that the dimensional gradient operator
$\nabla^{\prime}$ has the dimension of inverse of length):

$\displaystyle\mathcal{F^{\prime}}$
$\displaystyle=\frac{1}{V_{m}^{\prime}}\int_{\Gamma}\Big{[}f^{\prime}(c,\phi)+\kappa_{\phi}^{\prime}|\nabla^{\prime}\phi|^{2}+\kappa_{c}^{\prime}|\nabla
c|^{2}\Big{]}d\Gamma^{\prime}$ $\displaystyle\Rightarrow E_{c}\mathcal{F}$
$\displaystyle=\frac{1}{L_{c}^{3}V_{m}}\int_{\Gamma}\Big{[}E_{c}f(c,\phi)+(1/L_{c}^{2})\kappa_{\phi}^{\prime}|\nabla\phi|^{2}$
$\displaystyle+(1/L_{c}^{2})\kappa_{c}^{\prime}|\nabla
c|^{2}\Big{]}L_{c}^{3}d\Gamma,$ $\displaystyle\Rightarrow\mathcal{F}$
$\displaystyle=\frac{1}{V_{m}}\int_{\Gamma}\Big{[}f(c,\phi)+\frac{\kappa_{\phi}^{\prime}}{E_{c}L_{c}^{2}}|\nabla\phi|^{2}+\frac{\kappa_{c}^{\prime}}{E_{c}L_{c}^{2}}|\nabla
c|^{2}\Big{]}d\Gamma$ (15)

Thus we see that we get back the original form of Eq. (1) in terms of non-
dimensional variables when $\kappa_{c}^{\prime}$ and $\kappa_{\phi}^{\prime}$
are scaled by choosing a reference value of $\kappa=E_{c}L_{c}^{2}$.

The Cahn-Hilliard equation (Eq. 3) in dimensional form is converted to its
non-dimensional form as follows:

$\displaystyle\frac{\partial c}{\partial t^{\prime}}$
$\displaystyle=\nabla^{\prime}\cdot
M^{\prime}\nabla^{\prime}\Big{(}\frac{\delta\mathcal{F}^{\prime}}{\delta
c}\Big{)}$ $\displaystyle\Rightarrow\frac{1}{\tau_{c}}\frac{\partial
c}{\partial{t}}$
$\displaystyle=\frac{1}{L_{c}}\nabla\cdot\frac{1}{L_{c}}M^{\prime}\nabla\Big{[}\frac{\delta({E_{c}}\mathcal{F})}{\delta
c}\Big{]}$ $\displaystyle\Rightarrow\frac{\partial c}{\partial{t}}$
$\displaystyle=\nabla\cdot\frac{E_{c}\tau_{c}}{L_{c}^{2}}M^{\prime}\nabla\Big{(}\frac{\delta\mathcal{F}}{\delta
c}\Big{)}$ (16)

Thus we get back the non-dimensional Eq. (3) by scaling the dimensional
mobility $M^{\prime}$ with its characteristic value $L_{c}^{2}/E_{c}\tau_{c}$.
Using these conversion expressions and choosing appropriate values for the
reference variables $L_{c}$, $E_{c}$ and $\tau_{c}$, the dimensionless model
parameters used in the simulations can now be connected to their dimensional
counterparts [10].

Choosing $L_{c}=(V_{m}/N_{0})^{1/3}$ where $N_{0}$ is the Avogadro’s number
gives a reference length of $\sim$ $0.25\text{\,}\mathrm{nm}$. If reference
energy is taken as $k_{B}T_{c}$, it yields a value of $E_{c}=1.863\times
10^{-20}~{}$\mathrm{J}$=$0.116\text{\,}\mathrm{eV}$$ for
$T_{c}=$$1350\text{\,}\mathrm{K}$. Finally, a reference time can be obtained,
for example, by using a typical mobility value [22] of
$M^{\prime}\approx\dfrac{1.8\times 10^{-18}}{k_{B}T}$
$\text{\,}\mathrm{m}^{2}\mathrm{/}\mathrm{J}\mathrm{s}$, which for $T=$
$800\text{\,}\mathrm{K}$ yields $\tau_{c}=L_{c}^{2}/E_{c}M^{\prime}\approx 21$
ms. Dimensional values of relevant model parameters, along with their
conversion factors, are listed in Table LABEL:SimulationParameters, while the
corresponding values of surface and interfacial energies resulting from
different parameter sets are provided in Table LABEL:tab:IEdimensional.

Table 2: Simulation parameters (all energies are in _per atom_ basis).
Conversion factors from non-dimensional to dimensional form are based on
characteristic length $L_{c}=$0.25\text{\,}\mathrm{nm}$$, characteristic
energy $E_{c}=$0.116\text{\,}\mathrm{eV}$$ and characteristic time
$\tau_{c}=$21\text{\,}\mathrm{ms}$$.


††thanks: Communicating author

# Determination of the speed of light from the longitudinal modes beat
frequencies of an open-cavity helium-neon laser

Mohamad Taim<EMAIL_ADDRESS>Nicholas Kuder<EMAIL_ADDRESS>University of
Minnesota, Minneapolis, MN, USA

(August 27, 2024)

###### Abstract

We empirically determined the speed of light by measuring the variation in
longitudinal mode frequencies, or the beat frequencies, of an adjustable-
length, open-cavity helium-neon laser as a function of its cavity length. The
TEM00 mode lasing output of the laser was analyzed using a fast frequency
photodiode detector and a radio frequency spectrum analyzer. A Fabry-Perot
interferometer was used to monitor the intensity of the longitudinal modes and
we found that the phenomena of frequency pushing and pulling had little
effects on the beat frequency measurements. Plotting the reciprocal of the
beat frequency as a function of the change in cavity length, the speed of
light was found, by using linear weighted least squares regression, to be
$(2.997\pm 0.003)\times 10^{8}\,\mathrm{ms^{-1}}$. This value is $0.3\sigma$
away from the defined value of speed of light and is accurate to 1 part in
3200.

##  I Introduction

The speed of light, $c$, is an important fundamental constant in many
scientific fields, from electrodynamics to general relativity. $c$ is not just
about describing the propagation property of electromagnetic waves, but it
also happens to be the upper bound to the propagation speed of signals and to
the speeds of all material particles.

Large numbers of determinations have been made since Galilei Galileo first
attempted to prove the quantifiability of speed of light by using two lanterns
placed on top of two separate hills and measuring the time taken for the light
signals to cross the valley [1]. Most contemporary attempts on measuring the
speed of light employed variations of the two common experimental methods:
time of flight, and measurement of the frequency and wavelength of a radiation
[2]. Earlier papers described techniques like Foucalt’s rotating mirror method
[3], Kerr electro-optical shutter [4, 5], and microwaves [6].

Our experiment utilized the resonance condition of an open cavity, first done
by Essen and Gordon-Smith in 1947 with a cylindrical copper cavity resonator
as described in Ref. [7]. Their value of $c$ was accurate to 3 parts in
100,000 In 2010, the determination in Ref. [8] utilized a Helium-Neon (He-Ne)
laser to make a determination that was accurate to within 1 in 4100. We also
used a He-Ne laser system to determine $c$ by studying the longitudinal mode
frequencies of the lasing beam and its relationship with the length of the
laser cavity.

##  II Theory

A He-Ne laser system consists of a cavity in between two highly reflective
mirrors, or a Fabry-Perot cavity. A 99.9% reflective mirror is on one end and
a 99% reflective output coupler is on the other with the output coupler
letting through a little amount of light that becomes the laser beam. The
cavity is filled with a mixture of helium and neon and together, the gas-
filled cavity is called the medium.

When a DC current is pumped into the medium, the electrons in it excite into
higher levels of energy and the medium is said to undergo population
inversion, a process in which the number of excited atoms is more than the
number of ground atoms. These excited electrons then decay to lower levels of
energy and release photons. This is the spontaneous emission and the photons
are emitted in random directions and phases. Those that are parallel to the
horizontal axis of the laser tube get reflected by the two mirrors many, many
times. This leads to stimulated emission when the photons interact with other
excited electrons in the medium. The result is an avalanche of photons with
similar orientation, wavelength, and phase that comes out of the output
coupler as a bright, nearly monochromatic, single beam of laser.

The Fabry-perot cavity also imposes a resonance condition on the system in
that the laser will not operate with just any wavelength. The lasing output
will only have wavelengths, $\lambda$, that satisfy the equation

$\lambda=\frac{2nL}{k}\text{,}$ (1)

or, converted into mode frequency, $f$, in the equation

$f=\frac{kc}{2nL}\text{,}$ (2)

where $n$ is the index of refraction of the medium in which the light wave
travels, $L$ is the length of the cavity, and $k$ is an integer representing
the mode number [8].

As shown in Fig. 1, a laser beam has a Gaussian-shaped output [8]. This is due
to Doppler broadening, a spectrum of Doppler shifting effects caused by the
random directions of the photons in the medium.

Figure 1: The modes of a laser output under the medium gain curve.

Combined with the fact that the medium has a gain threshold, we can only see
two or three mode frequencies at a given time, and the difference between any
two of these mode frequencies is known as the beat frequency, $f_{b}$, given
as

$f_{b}=\frac{c}{2nL}\text{.}$ (3)

From this equation, assuming a fixed value for $n$, we can know $c$ accurately
by measuring and plotting only the two variables, $f_{b}$ and $L$.

##  III Experimental Setup & Methodology

###  III.1 Laboratory Setup

A He-Ne laser tube is mounted on an optical post and a $99.9\,\mathrm{\%}$
reflective output coupler (ThorLabs) with a radius of $45\,\mathrm{cm}$ is
secured on a translation stage (Micro-Controle), which is controlled by a
stepper motor (PhobOTronix), completing the open cavity of the laser system.
Between the output coupler and the laser tube is an adjustable iris to block
unwanted transverse modes.

The laser beam travels through a non-polarizing beam splitter, which directs
half of the beam to a fast frequency photodiode detector (Thorlabs model
PDA8A) and the other half into a scanning Fabry-Perot interferometer (SFPI)
(Tropel model 240). The signal from the photodetector is sent to a radio
frequency (RF) spectrum analyzer (Rigol model DSA815), which displays the beat
frequency. The beam from the SFPI is directed to another photodetector, which
is connected to an oscilloscope that will display the relative intensities of
the modes. Both the spectrum analyzer and the SFPI are assumed to be
calibrated correctly.

The measurements on beat frequency will be affected by two optical phenomena,
namely frequency pulling and pushing. Frequency pulling occurs when there are
varying index of refraction near the resonance transition. This causes the
difference in mode frequencies, depending on its position under the gain
profile curve, to have a smaller value [9]. In other words, the two mode
frequencies are being pulled toward the center. This effect is more pronounced
when we have mode frequencies that are not at similar intensity.

Meanwhile, frequency pushing is when we have an increase in gain, or increase
in refractive index, due to reasons that are myriad and complex, such as the
ambient temperature and pressure. This phenomenon increases the beat frequency
and also the total intensity [9]. We can prevent this by taking data only when
the total intensity is not changing.

The SFPI serves to monitor the relative intensities of the longitudinal modes.
The iris in between the Brewster window and the output coupler is used to
adjust the total intensity of the laser beam. These measures are taken to
minimize the uncertainty in beat frequency measurements.

The schematic diagram of the experiment is shown below in Fig. 2 [8].

Figure 2: The schematic diagram of the experimental setup.

###  III.2 Calibration

The stepper motor, controlled via a LabVIEW (version 2018) program, is used to
move the translation stage by a known number of steps. After the translation
stage translated the distance, a digital caliper, fixed on the translation
stage using tapes to minimize human error, is used to measure the total
distance moved by the translation stage. This is the change in laser cavity
length, $\Delta L$. This process is repeated until enough data points have
been collected. The change in length per step is then determined by the slope
of the data using linear regression.

##  IV Results & Data Analysis

In order to calibrate the translation stage and obtain a relationship between
the number of steps, $n$, and the translated distance, $\Delta L$, which
corresponds to the difference in length of the laser cavity, we measured
$\Delta L$ for a known value of $n$. We found from the plot in Fig. 3, for
$n=10$,

$\Delta L=2.5\times 10^{-5}\,\mathrm{m}$ (4)

The uncertainty in $\Delta L$ is $\sigma_{\Delta L}=\pm 1.0\times 10^{-5}$,
which is dictated by the tolerance of the digital caliper.

Figure 3: The plot for the calibration of the translation stage.

Using the iris to ensure that the mode is in TEM00 mode (single circular
beam), we then measured the beat frequencies as a result of varying the length
of the laser cavity. The values for the reciprocal of beat frequency,
$1/f_{b}$, as a function of $\Delta L$ is plotted in Fig. 4.

Figure 4: Reciprocal beat frequency as a function of the change in the open-
cavity length.

For each $\Delta L$, the value for $1/f_{b}$ was averaged over 100 sweep
cycles to minimize the uncertainty on frequency reading by an order of
magnitude.

Linear regression using the weighted least squares method was performed and
the data plotted in Fig. 4 yielded the final result for the speed of light in
air to be

$c=(2.997\pm 0.003)\times 10^{8}\,\mathrm{ms^{-1}}$ (5)

Our measured value is accurate to within $0.3\sigma$ of the defined value
$c=299,792,458\,\mathrm{ms^{-1}}$ [10]. As can be seen in Fig. 5, the linear
fit residuals show no structure, which indicates good data. Specifically, the
reduced chi squared is $\chi^{2}_{\mu}=0.922$ and the P-value is $P=0.58$.

Figure 5: Residuals from the linear fit.

The uncertainty in the final result can be explained by the phenomena
frequency pushing and pulling. Although we found no correlation between the
two phenomena and the beat frequency. Specifically, we found that the relative
intensities of the two modes observed on the SFPI oscilloscope had little
affect on the value of beat frequency on the spectrum analyzer. Further
investigation also found little effect on beat frequency caused by the total
intensity, which we controlled by the iris. However, we managed to get good
data in the sense that we measured the beat frequency only when the total
intensity and the relative intensity of the two modes are similar from one
measurement to the next.

In future experiments, one can reduce the uncertainty caused by frequency
pulling by implementing an electronic feedback system that can make slight
adjustments to the cavity length in order to maintain similar relative
intensities at all time. Other sources of error like the misalignment of the
optical equipment, and the thermal expansion of the optical table can also be
looked into.

##  V Conclusions

by varying the length of the open cavity of a He-Ne laser, we determined the
speed of light to be $c=(2.997\pm 0.003)\times 10^{8}\,\mathrm{ms^{-1}}$. This
value is $0.3\sigma$ away from the defined value
$c=299,792,458\,\mathrm{ms^{-1}}$. We found that the phenomena frequency
pulling and pushing had little to no effects on our measurements and thus,
negligible uncertainties on our final value of $c$.

###### Acknowledgements.

Our thanks to Professor Kurt Wick for his guidance and enthusiasm throughout
the research; to Professors Dan Dahlberg and Daniel Cronin-Hennesy for their
lectures; to Kevin Booth for the constant help with various problems. We
express gratitute to the Department of Physics for the equipments and the
permission to run the experiment.

## References

  * Foschi and Leone [2009] R. Foschi and M. Leone, Measurement of the velocity of light, and the reaction times, Perception 38, 1251 (2009).
  * Bradley  _et al._ [1972] C. C. Bradley, G. J. Edwards, D. J. E. Knight, W. R. C. Rowley, and P. T. Woods, Towards a new determination of the speed of light, Phys. Bull. 23, 15 (1972).
  * Domkowski  _et al._ [1972] A. J. Domkowski, C. B. Richardson, and N. Rowbotham, Measurement of the speed of light, Am. J. Phys. 40, 910 (1972).
  * Palmer and Spratt [1954] C. H. Palmer and G. S. Spratt, A laboratory experiment on the velocity of light, Am. J. Phys. 22, 481 (1954).
  * Rogers  _et al._ [1969] J. Rogers, R. McMillan, R. Pickett, and R. Anderson, A determination of the speed of light by the phase-shift method, Am. J. Phys. 37, 816 (1969).
  * Bol [1950] K. Bol, A determination of the speed of light by the resonant cavity method, Phys. Rev. 80, 298 (1950).
  * Essen  _et al._ [1948] L. Essen, A. C. Gordon-Smith, and C. G. Darwin, The velocity of propagation of electromagnetic waves derived from the resonant frequencies of a cylindrical cavity resonator, Proc. R. Soc. A 194, 348 (1948).
  * D’Orazio  _et al._ [2010] D. J. D’Orazio, M. J. Pearson, D. S. J. T. Schultz, M. W. Best, K. M. Goodfellow, R. E. Scholten, and J. D. White, Measuring the speed of light using beating longitudinal modes in an open-cavity HeNe laser, Am. J. Phys. 78, 524 (2010).
  * Lindberg [1999] Å. M. Lindberg, Mode frequency pulling in He–Ne lasers, Am. J. Phys. 67, 350 (1999).
  * [10] E. Tiesinga, P. J. Mohr, D. B. Newell, and B. N. Taylor, The 2018 CODATA recommended values of the fundamental physical constants, https://physics.nist.gov/cgi-bin/cuu/Value?c.