Theory for the single-particle dynamics in glassy mixtures with particle size swaps
Grzegorz Szamel







Department of Chemistry,
Colorado State University, Fort Collins, CO 80523




We present a theory for the single-particle dynamics in binary mixtures 
with particle size swaps. The general structure of the theory follows that 
of the theory for the collective dynamics in binary mixtures with particle size swaps,
which we developed previously [G. Szamel, Phys. Rev. E 98, 050601(R) (2018)].
Particle size swaps open up an additional relaxation channel, which
speeds up both the collective dynamics and the single-particle dynamics. 
To make explicit predictions, we resort to a factorization approximation 
similar to that employed in the mode-coupling theory of glassy dynamics. 
We show that, like in the standard mode-coupling theory, the single-particle 
motion becomes arrested at the dynamic glass transition predicted by
the theory for the collective dynamics. We compare the non-ergodicity 
parameters predicted by our mode-coupling-like approach for the an equimolar
binary hard sphere mixture with particle size swaps with the non-ergodicity 
parameters predicted by the standard mode-coupling theory for the same system 
without swaps. Our theory predicts that the “cage size” is bigger in the
system with particle size swaps.


Theory for the single-particle dynamics in 
glassy mixtures with particle size swaps
    Grzegorz Szamel
    December 30, 2023
====================================================================================





§ INTRODUCTION


Computer simulations are an essential tool in the quest for the understanding of 
glassy dynamics and the glass transition <cit.>. 
Their usefulness comes from the fact that they can access particles' 
positions and thus are able to provide much
more detailed picture of the dynamics than even most ingenious experiments. 
However, for many years computer simulations had to cope with the “glass ceiling”
<cit.>: since timescales accessible in simulations are orders
of magnitude shorter than experimental timescales, computer simulations utilizing
well equilibrated systems were only feasible in the temperature/density region
where the dynamics are only 4-5 decades slower than in normal stable liquids. 

Recently, a creative extension of the particle size swaps Monte Carlo,
that was first used in the context of glassy systems' simulations almost 30 years
ago <cit.>, was introduced. We should recall that since the original study
of Ref. <cit.>, the particle size swaps Monte Carlo method has been tried 
a few times <cit.>, but with very limited success. Then, Berthier and 
collaborators <cit.> turned the problem on its head 
by introducing several new glassy 
systems for which the particle size swaps Monte Carlo method speeds up the 
equilibration by many orders of magnitude. It is now possible to equilibrate 
glassy systems at temperatures comparable to and even lower than the laboratory glass
transition temperature. This allowed for the investigation of 
static properties of glassy systems at realistic conditions 
<cit.>.

The success of the particle size swaps method opened two sets of theoretical questions: 
questions about the method itself and questions about our fundamental understanding 
of the glass transition. The latter questions are concerned with the explanation
of the success of the particle size swaps method. While it is intuitively plausible
that particle size swaps open a new relaxation channel and thus should speeding up
the dynamics, the amount of the speed-up is quite significant and calls for 
a more detailed explanation. Additionally, one would like to understand theoretically
when does the particle size swaps method work and when it does not. Finally, 
one would like to describe the eventual slowing down of the particle size swaps
Monte Carlo method and compare it with the slowing down of local dynamics, either
Monte Carlo, or Brownian, or Newtonian dynamics.   

Out of the second set of questions opened by the success of the particle size 
swaps Monte Carlo, probably the most crucial one is concerned with the importance  
of cooperative processes operating on length scales that increase upon approaching 
the glass transition compared to the importance of local dynamics' slowing down 
<cit.>. Additional questions are concerned with the relevance 
of configurational entropies calculated theoretically by introducing order parameters
that explicitly forbid or allow particle size swaps 
<cit.>.

In spite of these interesting questions raised by the success of the particle size 
swaps Monte Carlo method, the theoretical analysis of this method is in its infancy. 
Three approaches have been proposed. First, Ikeda et al. 
<cit.> proposed to use the close relationship between
the dynamics of mean-field models and the replica theory-based description of
these models. Specifically, they used the correspondence between the dynamic 
transition predicted by the dynamic (mode-coupling-like) theory and identified with the 
dynamic crossover observed in simulations and the dynamic transition predicted by 
the replica approach. They argued that in the presence of particles size swaps
one has to adopt a more general structure of the ansatz for the 
inter-replica correlation functions, which are the order parameters within the
replica approach. They showed that the new ansatz leads to a shift of the
dynamic transition towards larger volume fractions, in qualitative agreement
with the simulations.

Next, Brito et al. <cit.>, in the context of continuous
distribution of particles' diameters, proposed to treat the diameters as additional 
variables, endowed with their own equations of motion (stochastic in the case of 
the Monte Carlo dynamics and deterministic continuous time equations in the case of
the Newtonian-like dynamics). This allowed Brito et al. to put forward some 
general arguments, which imply that allowing particle size swaps should lead to the 
decrease of the onset temperature for the glassy dynamics. However, their arguments 
do not lead to any quantitative prediction for this change nor to any  
prediction for the actual dynamical behavior with particle size swaps. Nevertheless,
Brito et al.'s idea of treating particles' diameters as additional variables
is extremely useful and it has already lead to generalizations of the original
particle size swaps Monte Carlo method <cit.>. 

We developed a dynamic theory for the acceleration due to the particle size
swaps <cit.>. 
The general structure of our theory agrees with the physical intuition:
particle size swaps act as an additional relaxation channel, which speeds up
the dynamics. To make explicit predictions we had to approximate the so-called 
irreducible memory functions, which are the main unknown theoretical quantities
in our approach. To this end we used a factorization approximation of the type
utilized in the mode-coupling theory of glassy dynamics and the glass
transition <cit.>. We calculated an approximate glass transition
phase diagram  for an equimolar binary hard-sphere mixture. The presence of particle 
size swaps shifts the dynamic glass transition towards higher volume fractions.
The shift increases with increasing ratio of the hard-sphere diameters; 
it saturates at about 4% at the diameter ratio of approximately 1.2.
We would like to emphasize that in our theory the set of order parameters
does not depend on whether particle size swaps are allowed or not. It is the
set of self-consistent equations for the order parameters that changes. This
is in contrast with the approach of Ikeda et al. 

Here we extend the previously presented theory to the case of single-particle
motion. This extension will allow us to study the frozen-in part of the self-van 
Hove function
(albeit in the reciprocal space) and compare the localization lengths with
and without particle size swaps. Notably, the set of order parameters for the 
single-particle motion does depend on whether particle size swaps are allowed or not.
This is due to the fact that without particle size swaps there cannot be any
single-particle cross-correlations. 

Generally speaking, the relation of the structure of the theory for the single-particle
motion and the previously derived theory for the collective dynamics is very similar
to the relation between the same theories within the standard mode-coupling 
theory of the glass transition. In particular, in the present case involving
particle size swaps and in the standard theory describing local dynamics,
single-particle motion becomes arrested at the dynamic glass transition 
predicted by the theory for collective dynamics.

The paper is organized as follows. In the next section we discuss the model
system with particle size swaps that we analyze and introduce the self-intermediate
scattering functions that quantify the single-particle dynamics. In Sec. 3 we 
express these functions in terms of the so-called memory functions. In Sec. 4
we introduce the main approximation of the theory, the factorization approximation,
and use it to express the memory functions in terms of the intermediate
scattering functions. In Sec. 5 we present numerical results for the non-decaying
parts of the intermediate scattering functions in the glassy phase predicted by the 
theory. We end with some discussion in Sec. 6. 



§ MODEL: BINARY MIXTURE WITH PARTICLE SIZE SWAPS


We consider the model investigated in Ref. <cit.>, i.e. a binary 
mixture. This is the simplest model that allows one to investigate the influence 
of the particle size swaps on the dynamics. Two arguments can be made against
using a binary mixture. First, computer simulations studies <cit.> showed 
that particle size exchanges in systems with a continuous polydispersity result
in the largest speed-up of the dynamics. Second, a system with a continuous 
polydispersity allows one to treat the particles' diameters as additional variables 
that evolve according to their own equation of motion <cit.>.
The distribution of the particles diameters is then enforced by a diameter-dependent 
chemical potential-like term, which plays the role of an one-particle external 
potential for the diameter variable. The evolution of particles positions and 
sizes could then be implemented in a very similar way. However, the price that one has 
to pay for this uniformity of the evolution equations is that one has to deal with
a static problem with infinitely many components. This makes explicit calculations
a bit daunting.  

To investigate the single-particle dynamics we consider a binary mixture consisting 
of N particles in volume V. We assume that one of these particles, for convenience
the particle number 1, is tagged and focus our attention on the motion of this particle.
All the particles, including the tagged particle, can be of type A or B.
These particles types differ by size. Since any particle 
can change its type, the state of particle i is determined by its position,
𝐫_i, and type indicated by a binary variable σ_i,
with σ_i=1 corresponding to A and σ_i=-1 corresponding to B. 
The composition of the system is specified by the difference of the chemical
potentials of particles of type A and B, Δμ. We note that for 
a given chemical potential difference, the
composition depends also on the number density n=N/V and the temperature T
of the system. We will assume that these three parameters result in 
concentrations x_A and x_B=1-x_A, and that these concentration are constant while
the density or the temperature of the system varies. We note that this implies that
the tagged particle can be of type A with probability x_A and of type B with
probability x_B. In practical calculations 
we will restrict ourselves to equimolar mixtures (for which some formulas simplify).
We assume that the “normal” dynamics of the system is Brownian, i.e., that 
each particle moves under the combined influence of thermal noise and interparticle 
forces. The forces are derived from a spherically symmetric potential, which depends 
on the particle type, V_σ_i σ_j(r_ij), where 
r_ij=|𝐫_i-𝐫_j| is the distance between particles i and j.
In addition to the Brownian motion in space, the particles can change their type (size). 
We assume that each particle can change its type independently of the type changes
of other particles. This is analogous to the single-spin-flip dynamics of spin
systems. In contrast, in the initial computational 
studies <cit.> particle size swaps were done 
in such a way that the number of particles of each
size was conserved. This procedure is analogous to the so-called Kawasaki dynamics 
of spin systems. It is convenient in simulational studies because it allows one to 
maintain easily a specific composition of the system. Later investigations
<cit.> showed that replacing the Kawasaki-type procedure by 
single particle size swaps is a relatively mild change 
and that both procedures lead to very similar results. 

The above described model corresponds to the following equation of motion 
for the N-particle distribution,
P_N(𝐫_1,σ_1, ..., 𝐫_N,σ_N;t), abbreviated below as P_N,

    ∂_t P_N =
    Ω_sw
    P_N ≡(Ω + δΩ_sw) P_N.

Here the evolution operator Ω_sw consists of two parts
describing two relaxation channels, the part 
describing Brownian motion of particles,

    Ω = D_0 ∑_i ∂_𝐫_i·(∂_𝐫_i-β𝐅_i)

and the part describing particle size swaps,

    δΩ_sw = 
    - τ_sw^-1∑_i (1-S_i) w_i .

In Eq. (<ref>) D_0 is the diffusion coefficient of an isolated particle,
D_0=k_BT/ξ_0, with ξ_0 being the friction coefficient of an isolated particle,
β=1/k_B T, and 
𝐅_i is the total force on particle i, 
𝐅_i=∑_j≠ i𝐅_ij = - ∑_j≠ i∂_𝐫_i
V_σ_iσ_j(r_ij). In Eq. (<ref>)
τ_sw is the rate of attempted particle size swaps, S_i is the swap 
operator, S_i σ_i = -σ_i, and w_i is the factor ensuring that the 
detailed balance condition, (1-S_i)w_i P_N^eq = 0, is satisfied,
with P_N^eq being the equilibrium distribution,
P_N^eq∝exp(-β∑_i≠ j V_σ_iσ_j(r_ij) + 
∑_i β/2Δμσ_i ). 
The factor w_i depends on the way particle size swaps are attempted. In practical
applications one typically uses Metropolis criterion for accepting attempted swaps. 
It should be emphasized that while the interactions influencing particles' motion
in space are pairwise-additive, the factor w_i typically is not and it depends on
the whole neighborhood of particle i. 

The focus of the present investigation are the elements of the matrix
of the self-intermediate correlation functions, i.e. the 
single-particle (tagged particle) density correlation functions,

    F_αβ^s(q;t)  
    = <n_α^s(𝐪)e^Ω_swt 
    n_β^s(-𝐪)>,

where n_α^s is the Fourier transform of the microscopic 
tagged particle density, 

    n_α^s(𝐪) 
    =   
    1+σ_1(δ_α A-δ_α B)/2 
    e^-i𝐪·𝐫_1.


To develop a theory for the dynamics of F_αβ^s(q;t) we will also
need the collective density correlation functions,

    F_αβ(q;t)  
    = <n_α(𝐪)e^Ω_swt 
    n_β(-𝐪)>,

where n_α, α=A,B, 
are the Fourier transforms of the normalized microscopic densities
of particles of type α,

    n_α(𝐪) 
    =  1/√(N)∑_i>11+σ_i(δ_α A-δ_α B)/2 
    e^-i𝐪·𝐫_i.

In Eqs. (<ref>,<ref>) 
and in the following equations the standard conventions apply: 
<…> denotes the semi-grand canonical ensemble average
over P_N^eq, the equilibrium probability distribution
stands to the right of the quantity being averaged, and all operators
act on it as well as on everything else.

We should emphasize the contrast between the single-particle and collective density 
correlation functions. In the absence of particle size swaps the matrix
of the self-intermediate correlation functions is diagonal and the following
notation is used F_α^s(q;t)≡ F_αα^s(q;t), α=A,B.
In the presence of particle size swaps the tagged particle can change its type,
which leads to the appearance of the off-diagonal terms in the  matrix
of the self-intermediate correlation functions. In contrast, 
the set of collective density correlation functions is the same regardless of the
presence of particle size swaps. 

In the approach of Ikeda et al., the order parameters are the replica theory
analogues of the single-particle density correlations (physically, the long-time
limits of the time-dependent single-particle density correlation functions should
be equal to the replica theory single-particle correlations; note that a truly
consistent dynamic and static description of the glass transition is still
lacking, except in infinite spatial dimension). Thus, one could argue that it
is natural that in that approach a different set of order parameters is used
depending on the presence of particle size swaps. 



§ MEMORY FUNCTION REPRESENTATION



To derive the memory function representation of the tagged particle
density correlations we use the standard projection operator procedure
<cit.>. First, we define a projection operator on the tagged particle
density subspace, 𝒫^s,

    𝒫^s = …∑_α. n_α^s(-𝐤)>x_α^-1< n_α^s(𝐤) . … .

and the orthogonal projection, 𝒬,

    𝒬^s = ℐ - 𝒫^s ≡ℐ - …∑_α. n_α^s(-𝐤)>x_α^-1< n_α^s(𝐤) . … .

Note that <n_α^s(𝐤) n_β^s(-𝐤)> ≡ 
F_αβ^s(𝐤;t=0)=x_αδ_αβ.

Next, we use the standard projection operator identity <cit.>

    1/z-Ω_sw = 1/z - Ω_sw𝒬 +
    1/z - Ω_sw𝒬Ω_sw𝒫1/z-Ω_sw.

to express the Laplace transform of the time derivative of the tagged particle 
density correlation function,
zF_αβ^s(q;z)-F_αβ^s(q;t=0), 
in terms of the tagged particle reducible memory function,

    zF_αβ^s(q;z)-F_αβ^s(q;t=0)
    = <n_α^s(𝐪) Ω_sw1/z-Ω_sw n_β^s(-𝐪) > 
    
       =   <n_α^s(𝐪) Ω_sw( 𝒫^s +  𝒬^s )
    1/z-Ω_sw n_β^s(-𝐪) >
    
       =   ∑_γ<n_α^s(𝐪) Ω_sw n_γ^s(-𝐪)>x_γ^-1< n_γ^s(𝐪) 1/z-Ω_sw 
    n_β^s(-𝐪) > 
    
           + 
    ∑_γ<n_α^s(𝐪) Ω_sw𝒬^s  
    1/z - Ω_sw𝒬^s𝒬^s Ω_sw n_γ^s(-𝐪)>x_γ^-1< n_γ^s(𝐪)1/z-Ω_sw
    n_β^s(-𝐪) >.

In the second line of Eq. (<ref>) we identify the tagged particle 
frequency matrix,

    O_αβ^s(q)    =    
    - <n_α^s(𝐪) Ω_sw n_β^s(-𝐪)> =
    D_0 𝐪^2 x_αδ_αβ +
    τ_sw^-1<n_α^s(𝐪)
    (1-S_1) w_1 n_β^s(-𝐪)>
    
       =   
    D_0 𝐪^2 x_αδ_αβ +  
    1/2τ_sw< w_1 > (2δ_αβ-1),

and in the third line we identify the matrix of the tagged particle reducible  
memory functions,

    M_αβ^s red(q;z)    =   <n_α^s(𝐪) Ω_sw𝒬^s 
    (z-𝒬^sΩ_sw𝒬^s)^-1
    𝒬^s Ω_sw n_β^s(-𝐪)>.


We note that, like in the collective dynamics problem <cit.>, 
the tagged particle frequency matrix can be decomposed into the part corresponding
to “normal” (i.e. Brownian) and the part describing particle size swaps.
This separation can be conveniently expressed by introducing 3-dimensional vectors 
𝗏_α, α=A,B, and a 3x3 matrix 𝖮^s,

    O_αβ^s(q)    =   𝗏_α^T𝖮^s(q) 𝗏_β

Here 𝗏_A^T=(1,0,1), 𝗏_B^T=(0,1,-1), 
𝖮_11^s=D_0q^2 x_A, 𝖮_22^s=D_0q^2 x_B, 
𝖮_33^s=(1/2τ_sw)<w_1> and 
𝖮_ab=0 for a≠ b. Thus, Brownian dynamics is represented by 
elements of matrix 𝖮^s with a,b ≤ 2 and the particle size swaps 
are represented by the 33 element of 𝖮^s. 

Next, we note the structure of the “vertexes” in the 
matrix of the tagged particle reducible memory functions,

    . 𝒬^s Ω_sw n_α^s(-𝐪)>    =   . D_0
    𝒬^s i 𝐪·( i𝐪 + β𝐅_1)
    1+σ_1(δ_α A-δ_α B)/2
    e^i𝐪·𝐫_1>
    
          
    - . 1/τ_sw𝒬^s 
    w_1 σ_1 e^i𝐪·𝐫_1>(δ_α A-δ_α B)
    
       =   ∑_a . π_a^s(-𝐪)> 𝗏_α a,

where functions π_a^s are defined as follows,

    π_1,2^s(𝐪)    =    
    D_0
    𝒬^s i𝐪·( i𝐪 - β𝐅_1)
    1±σ_1/2e^-i𝐪·𝐫_1,
    
    π_3^s(𝐪)    =    - 1/τ_sw𝒬^s 
    w_1 σ_1 e^-i𝐪·𝐫_1.

and 𝗏_α a is the ath element of vector 𝗏_α.

The structure of the vertexes uncovered in Eq. (<ref>) 
allows us to express the matrix of the tagged particle reducible
memory functions M^s red in terms of  3-dimensional vectors 
𝗏_α, α=A,B, 
and a 3x3 matrix 𝖬^s red,

    M_αβ^s red(q;z) 
       =   𝗏_α^T𝖬^s red(q;z) 
    𝗏_β,

where the matrix elements of 𝖬^s red read 

    𝖬_ab^s red(q;z) = 
    <π_a^s(𝐪) 
    (z-𝒬^sΩ_sw𝒬^s)^-1π_b^s(-𝐪)>.


Again, this representation allows us  
to separate “normal” (i.e. Brownian) dynamics represented by 
elements of matrix 𝖬^s red with a,b ≤ 2 
and the particle swaps represented by the 33 element
of 𝖬^s red. We note, however, that while matrix 𝖮^s
is diagonal, in general, matrix 𝖬^s red is not. Thus, 
in general, there will be terms describing the time-delayed couplings between 
the two relaxation channels. 
These couplings are described by elements 𝖬^s red with
a ≤ 2 and b=3, and a=3 and b≤ 2. 

Finally, we introduce the tagged particle irreducible evolution operator. This 
will allow us to introduce a 3x3 matrix 𝖬^s irr, 
whose elements are functions evolving with the so-called irreducible evolution operator.
For systems evolving with Brownian dynamics this step was introduced and justified 
by Cichocki and Hess <cit.>. Later, Kawasaki <cit.> argued that an
analogous construction should also be applied for a more general class
of systems evolving with stochastic dynamics, including spin systems evolving with 
spin-flip dynamics. We recall that our dynamics is a combination of Brownian
dynamics and particle size swaps, and the latter events are technically implemented 
in terms of single spin flips.

We follow the definition of the collective dynamics irreducible evolution 
operator <cit.> and define the tagged particle irreducible 
evolution operator as follows,

    Ω_sw^s irr   =   𝒬^sΩ_sw𝒬^s -
    ∑_αβ. 𝒬^s Ω n_α^s(-𝐪)>
    <n_α^s(𝐪)Ω n_β^s(-𝐪)>^-1<n_β^s(𝐪)Ω𝒬^s. 
    
           
    -
    . 𝒬^s δΩ_swσ^s(-𝐪)>
    <σ^s(𝐪)δΩ_swσ^s(-𝐪)>^-1<σ^s(𝐪)δΩ_sw𝒬^s. ,

where σ^s(𝐪)=σ_1e^-i𝐪·𝐫_1. 
With the help of matrix 𝖮^s Eq. (<ref>) 
can be re-written as

    Ω_sw^s irr   =   𝒬^sΩ_sw𝒬^s +
    ∑_a,b^2 . π_a^s(-𝐪)> (𝖮^s)^-1_ab<π_b^s(𝐪). +
    . 2 π_3^s(-𝐪)>(4𝖮^s)^-1_33<2π_3^s(𝐪). 
    
       =   𝒬^sΩ_sw𝒬^s +
    ∑_a,b. π_a^s(-𝐪)> (𝖮^s)^-1_ab(q)
    <π_b^s(𝐪). .

We note that, like in the case of the collective dynamics irreducible operator,
one-particle reducible parts of 𝒬^sΩ_sw𝒬^s are removed
separately for Brownian dynamics and particle size swaps. 

Next, we use
a formula analogous to Eq. (<ref>),

    1/z-𝒬^sΩ_sw𝒬^s   =   1/z - Ω_sw^s irr
           -  
    1/z - Ω_sw^s irr∑_a,b. π_a^s(-𝐪)> (𝖮^s)^-1_ab(q)
    <π_b^s(𝐪). 
    1/z-𝒬^sΩ_sw𝒬^s

and we express 𝖬^s red in
terms of matrix 𝖬^s irr whose elements are functions 
evolving with the irreducible evolution operator,

    𝖬^s red(q;z) =
    𝖬^s irr(q;z) - 
    𝖬^s irr(q;z) (𝖮^s)^-1(q)
    𝖬^s red(q;z),

where 

    𝖬_ab^s irr(q;z) = 
    <π_a^s(𝐪) (z-Ω_sw^s irr)^-1π_b^s(-𝐪)>.


The combination of Eqs. (<ref>-<ref>), (<ref>) 
and (<ref>-<ref>) constitutes the memory function representation
for the tagged particle density correlation functions. 
We note that up to now, our approach was formally exact. On the other hand, 
it could be said that we have just re-wrote the problem of calculating
the tagged particle density correlation functions as a problem of
calculating the elements of matrix 𝖬^s irr. Thus, the 
latter functions have become the main unknown objects of our theory.



§ FACTORIZATION APPROXIMATION


To proceed, we need to introduce some approximations. Here, as in the earlier
study of the collective dynamics, we follow the spirit of the mode-coupling theory 
for glassy dynamics and the glass transition <cit.> and 
use a factorization approximation. 
More precisely, we use a sequence of three approximations 
<cit.>.  
In the present case these approximations need to be slightly modified to 
take into account the focus on the tagged particle dynamics.  

First, in expression (<ref>) for the 
matrix elements of 𝖬^s irr 
we project functions π_a^s onto the subspace spanned by the part of the product 
of the tagged particle density and the density of the other particles
orthogonal to the tagged particle density. The resulting expression reads,

    𝖬^a irr_ab(𝐪; t) ≈∑_α,...,θ∑_𝐤_1,...,𝐤_x<π_a^s(𝐪) 
    n_α^s(-𝐤_1)n_β(-𝐤_2)> 
    
          ×<n_α^s(𝐤_1)n_β(𝐤_2)𝒬^s
    n_γ^s(-𝐤_3)n_δ(-𝐤_4)> ^-1
          ×<n_γ^s(𝐤_3)n_δ(𝐤_4)
    exp(Ω_sw^s irrt)
    n_ϵ^s(-𝐤_5)n_ζ(-𝐤_6)>
    
          ×<n_ϵ^s(𝐤_5)n_ζ(𝐤_6)𝒬^s
    n_η^s(-𝐤_7)n_θ(-𝐤_8)>^-1< n_η^s(𝐤_7)n_θ(𝐤_8)π_b^s(-𝐪)
    >.

In the case of Brownian dynamics with pairwise additive interactions the
first approximation is exact. In the present case of Brownian dynamics
with particle size swaps, Eq. (<ref>) constitutes an approximation, due
to the fact that factors w_i are typically not pairwise additive.  

Second, we factorize four-point dynamic 
correlation functions while replacing the tagged particle irreducible evolution 
operator by the original un-projected evolution operator, 

    <n_γ^s(𝐤_3)n_δ(𝐤_4)𝒬^s 
    exp(Ω_sw^s irr t) 
    𝒬^s n_ϵ^s(-𝐤_5)n_ζ(-𝐤_6)>
    
       ≈   <n_γ^s(𝐤_3) exp(Ω_sw t) 
    n_ϵ^s(-𝐤_5) >
    <n_δ(𝐤_4) exp(Ω_sw t) n_ζ(-𝐤_6) >
    
       ≡   <n_γ^s(𝐤_3;t) n_ϵ^s(-𝐤_5) >
    <n_δ(𝐤_4;t) n_ζ(-𝐤_6) >.

We note that Eq. (<ref>) is the crucial approximation of our approach.

Third, we use some additional approximations for static correlation functions. 
We factorize normalization factors
<n_α^s(𝐤_1)n_β(𝐤_2)𝒬^s
n_γ^s(-𝐤_3)n_δ(-𝐤_4)> ^-1,

    <n_α^s(𝐤_1)n_β(𝐤_2)𝒬^s
    n_γ^s(-𝐤_3)n_δ(-𝐤_4)>^-1
          ≈< n_α^s(𝐤_1)n_γ^s(-𝐤_3) >^-1< n_β(𝐤_2)n_δ(-𝐤_4) >^-1,

and we note that < n_α^s(𝐤_1)n_γ^s(-𝐤_3) >^-1
=x_α^-1δ_𝐤_1𝐤_3δ_αγ and
< n_β(𝐤_2)n_δ(-𝐤_4) >^-1 = 
δ_𝐤_2𝐤_4S_βδ^-1(k_2), with S_βδ(k)
being the partial static structure factor. We note that the vertices
originating from Brownian dynamics,
< n_η^s(𝐤_7)n_θ(𝐤_8)π_b^s(-𝐪)>,
b=1,2, can be calculated exactly, 

    < n_η^s(𝐤_7)n_θ(𝐤_8)π_b^s(-𝐪)
    > = nx_b/√(N)δ_𝐤_7+𝐤_8,𝐪𝐤_8·𝐪δ_η b∑_γ S_θγ(k_8)c_γ b(k_8)
          b=1,2.

The calculation of the vertex originating from particle size swaps,
< n_η^s(𝐤_7)n_θ(𝐤_8)π_3^s(-𝐪)>,
is a bit more involved. Using the mixture version of the convolution 
approximation we initially obtain the following expression

    < n_η^s(𝐤_7)n_θ(𝐤_8)π_3^s(-𝐪)> 
    = -(2τ_sw√(N))^-1δ_𝐤_7+𝐤_8,𝐪(δ_η A-δ_η B)
    
          ×[R_θ(k_8)-1/2(Q(δ_θ A-δ_θ B)
    +<w_1>)
    -< w_1 >x_η^-1(S_ηθ(k_8)-x_θδ_ηθ)]

where 

    R_θ(q) = < n_θ(𝐪)
    w_1 e^i𝐪·𝐫_1> 
        and     Q = < σ_1 w_1 >

To simplify expression (<ref>) we factor out 
<w_1> from R_θ and Q,

    R_θ(q) ≈<w_1 >
    ∑_λ S_θλ(q)  
        and     Q ≈< w_1 > < σ_1 >
    = < w_1 > (x_A-x_B)

Using Eq. (<ref>) in Eq. (<ref>) we obtain 

    < n_η^s(𝐤_7)n_θ(𝐤_8)π_3^s(-𝐪)> 
    = -(2τ_sw√(N))^-1δ_𝐤_7+𝐤_8,𝐪(δ_η A-δ_η B)
    
          ×<w_1> [∑_λ S_θλ(k_8)
    -x_Aδ_θ A - x_Bδ_θ B
    - x_η^-1(S_ηθ(k_8)-x_θδ_ηθ)].


The final result is a set of approximate expressions for the
elements of matrix 𝖬^s irr in terms of integrals of products
of the tagged particle and collective density correlation functions. Here we list 
these expressions for the equimolar mixture, x_A=x_B=0.5, which was considered
in Ref. <cit.>, in the thermodynamic limit, 
N→∞, V→∞, N/V=n=const. 
The expressions for 𝖬^s irr_ab,
a,b=1,2 are the same as those derived in the standard mode-coupling theory,

    𝖬^sirr_ab(𝐪; t)    ≈    nD_0^2 
    ∫d𝐤/(2π)^3∑_βδ𝐪·𝐤
    c_aβ(k) F^s_ab(|𝐪-𝐤|;t)
    
          ×
    F_βδ(k;t)c_δ b(k)
    𝐪·𝐤.

In Eq. (<ref>) and in the following equations, 
c_α 1≡ c_α A and c_α 2≡ c_α B.
The remaining matrix elements originate from particle size swaps,

    𝖬^s irr_a3(𝐪; t)    ≈    
    -nD_0 < w_1 >/2τ_sw∫d𝐤/(2π)^3∑_βγδ𝐪·𝐤
    c_aβ(k) F^s_aγ(|𝐪-𝐤|;t)
    
          ×
    F_βδ(k;t)(c_δ A(k)-c_δ B(k)) = 
    𝖬^s irr_3a(𝐪; t)       a=1,2,


    𝖬^s irr_33(𝐪; t)    ≈    
    n< w_1 >^2/(2τ_sw)^2∫d𝐤/(2π)^3∑_αβγδ(c_Aβ(k)-c_Bβ(k))
    F^s_αγ(|𝐪-𝐤|;t)
    
          ×
    F_βδ(k;t)(c_δ A(k)-c_δ B(k)).


The exact memory function representation defined through Eqs. 
(<ref>-<ref>), (<ref>) and (<ref>-<ref>), 
combined with approximate form of the elements of matrix 𝖬^s irr 
in terms of integrals of products of the tagged particle and collective density 
correlation functions, Eqs. (<ref>-<ref>), constitute our 
mode-coupling theory for the single-particle motion. Numerical solution of these
time-dependent equations looks quite a bit more difficult than the solution
of the standard mode-coupling equations for mixtures. The reason for this supposition
is that in the present case the relation between the memory function and the
correlation functions is more complicated.

We note that, as in the standard mode-coupling theory, the necessary inputs into
the theory for the single-particle dynamics are collective intermediate scattering
functions, which in the present case should be obtained from the theory for
the collective dynamics with particle size swaps developed in Ref. <cit.>.



§ NUMERICAL RESULTS FOR THE ARRESTED PHASE


First, let us recall the analysis 
presented in Ref. <cit.>. We started from equations for the 
time-dependent collective density correlation functions derived 
and we assumed that as the density increases and/or 
the temperature decreases, collective density correlations develop plateaus, 
and that at the transition nonzero long-time limits F_αβ(q; ∞) appear 
discontinuously. These assumptions allowed us to derive a set of self-consistent 
equations for F_αβ(q; ∞). We analyzed these self-consistent equations 
numerically for an equimolar binary hard-sphere mixture using as the static input 
equilibrium correlation functions obtained from the Percus-Yevick (PY) closure 
<cit.>. We obtained a fluid-glass phase diagram. 
We found that for a system with particle size swaps the volume fraction 
at the dynamic glass transition increases with increasing particle diameter ratio 
up to the ratio of about 1.2, where the relative increase is about 4%, and then
saturates. In agreement with Ref. <cit.>, the location of the 
dynamic glass transition for an equimolar mixture without particle size swaps 
depends very weakly on the ratio of particle diameters for the ratio smaller than 1.5.

Here, we focus on the non-decaying parts of the tagged particle density
correlations F_αβ^s(q; ∞). Their appearance signals the arrest of
the single-particle motion. We recall that in the standard mode-coupling theory
the dynamic glass transition, i.e. the appearance of the frozen-in
collective density fluctuations, coincides with the arrest of the single-particle
motion. We note that this prediction is non-trivial. Generally speaking, if one
considers single-particle motion in a fluid with some frozen-in disorder, the
arrest of the single-particle motion is only predicted if the disorder is strong enough.
It follows from an explicit calculation that the frozen-in density fluctuations 
predicted by the standard mode-coupling theory at the dynamic glass transition 
are strong enough to arrest the single-particle motion (under the assumption that 
the latter process is described by the mode-coupling theory for the single-particle 
motion). 

First, we derive the self-consistent equations for the non-decaying parts of the 
tagged particle density correlations F_αβ^s(q; ∞). We start 
from the combination of Eqs. (<ref>-<ref>), (<ref>), 
(<ref>-<ref>) and (<ref>-<ref>)
and obtain

    F_αβ^s(q;∞) - x_αδ_αβ(q) = 
           
    - ∑_γδ𝗏^T_α𝖮^s(q) 
    [𝖬^s irr(q,∞)]^-1𝖮^s(q) 𝗏_γ
    x_γ^-1δ_γδ F_δβ^s(q;∞),

where 𝖬^s irr(q,∞) is given by 
Eqs. (<ref>-<ref>) with non-zero long-time limits 
F_αβ(q;∞) and F_αβ^s(q;∞) 
substituted at their right-hand-sides.

We solved self-consistent equations (<ref>) numerically <cit.>, 
using as the static 
input equilibrium correlation functions obtained from the Percus-Yevick (PY) closure 
<cit.> and using previously obtained non-decaying parts
of the collective density fluctuations F_αβ(q;∞). We found that,
as is the case in the standard mode-coupling theory, the single-particle
motion becomes arrested at the location of the dynamic glass transition. 

In the first two sets of figures below, Figs. <ref>-<ref>, 
we compare the results obtained from the present theory for 
an equimolar binary hard-sphere mixture with particle size swaps with the 
results of the standard mode-coupling theory for the same system with Brownian dynamics
only. In both cases we present results at the dynamic glass transition, respectively 
the transition predicted by the present theory and that predicted by the standard
mode-coupling theory. In all the figures we present results for the 
diameter ratio at which the difference between these transitions is approximately
the largest, d_L/d_S=1.2.

In Figs. <ref> we compare the non-decaying parts of the 
collective density correlations for the larger particles, F_AA(q; ∞),
with the corresponding partial structure factors S_AA(q) at the same
volume fraction. We see that in both cases the wavevector dependence of 
F_AA(q; ∞) follows that of S_AA(q). The main difference between 
the two results is seen in the zero wavevector limit. 
The present theory predicts lim_q→ 0 F_AA(q; ∞) is significantly smaller 
than lim_q→ 0 S_AA(q), in contrast with the 
standard mode-coupling binary mixture result, 
lim_q→ 0 F_AA(q; ∞) ≈lim_q→ 0 S_AA(q). 
This qualitative difference originates from the particle size swaps.



In the next set of figures we compare the non-decaying parts of the reduced, diagonal 
collective density correlations F_αα(q; ∞)/S_αα(q), 
α=A,B and tagged particle density fluctuations, F_αα^s(q; ∞),
obtained from the present theory, Figs. <ref>a-b, and from the
standard mode-coupling theory, Figs. <ref>a-b. Now we clearly see the 
difference of the small wavevector behavior of the frozen-in collective density
correlations discussed above. In addition, we see the qualitative difference 
between the non-decaying parts of the diagonal tagged particle 
density correlations F_αα^s(q; ∞). The present theory predicts
lim_q→ 0 F_αα^s(q; ∞)/x_α < 1 whereas within the standard 
mode-coupling approach the conservation of the number of the particles of each 
type leads to lim_q→ 0 F_α^s(q; ∞)/x_α = 1.






Following the spirit of Ikeda et al. <cit.> we can
use the small wavevector limits of the frozen-in tagged particle 
density correlations F_αα^s(q; ∞) to define an order parameter
for particle size swaps in the glass,

    q= lim_q→ 0[F_AA^s(q; ∞)-F_AB^s(q; ∞)-F_BA^s(q; ∞)+F_BB^s(q; ∞)
    ],

where use the same symbol q as in Refs. <cit.>
but caution the reader no to confuse this quantity with the wavevector q.
The zero value of the order parameter, q=0 indicates the complete lack 
of particle type correlations in the glass whereas the maximum value, q=1,
indicates the absence of the particle size swaps in the glass.

In Fig. <ref>a we show the volume fraction dependence of the 
particle size swaps order parameter for the 
diameter ratio at which the difference between these transitions is approximately
the largest, d_L/d_S=1.2. On physical grounds one would expect the particle size 
swaps to become less probable with increasing volume fraction. The present theory
confirms this expectation: we observe q increasing monotonically starting 
from a nonzero value q≈ 0.376 at the dynamic glass transition.
Qualitatively, the volume fraction dependence of the particle size swaps order parameter
predicted by the present theory agrees with the prediction of the replica approach
of Ikeda and Zamponi <cit.>. 



Within the standard mode-coupling theory, 
the small wavevector behavior of the frozen-in tagged particle 
density correlations F_α^s(q; ∞) allows one to define the 
long-time limit of the mean square displacement in the glass, which can
be interpreted as a quantification of the “cage size”,

    F_α^s(q;∞) = x_α( 1 - 1/2A_α q^2/2
    + ... )   where   A_α=
    lim_t→∞<|𝐫_1^α(t)-𝐫_1^α(0)|^2>.
 

In the present case of a system with particle size swaps we generalize this 
definition as

    F_αβ^s(q;∞) = √(x_α x_β)( 1 - 1/2A_αβ q^2/2
    + ... ) 
    
                     where    A_αβ=
    lim_t→∞<|𝐫_1^α(t)-𝐫_1^β(0)|^2>.
 

In Fig. <ref>b we show the volume fraction dependence of the
A_αβ and A_α for the 
diameter ratio at which the difference between these transitions is approximately
the largest, d_L/d_S=1.2. As expected on physical grounds, all the quantities 
characterizing “cage size” decrease with increasing volume fraction. For 
a given volume fraction, the “cage size” for the system with particle size swaps is
larger than that for the system without swaps. This is consistent 
with simulational results of Ninarello et al <cit.>. 
Furthermore, as expected on physical grounds, the cage size depends on the 
particle type and is smaller for the larger particles. 
Again, qualitatively, the volume fraction dependence of the “cage size” 
predicted by the present theory agrees with the prediction of the replica approach
of Ikeda and Zamponi <cit.>. However, the important difference is
that the latter theory uses only a single cage size parameter for both types of the
particles, both in the 
system with particle size swaps and in the system without swaps.




§ DISCUSSION


We presented here a theory for the single-particle motion in binary
mixtures with particle size swaps. The present theory is an extension of the
previously developed theory for the collective dynamics and the dynamic glass
transition in binary mixtures with particle size swaps. Like the latter approach,
the present theory is based on combination of a formally exact memory function
representation of the tagged particle density correlation function and a 
factorization approximation for a four-point correlation function evolving with 
the irreducible dynamics. Thus, it belongs to the class of mode-coupling-like
theories. 

The present theory predicts that the single-particle motion becomes localized
at the dynamic glass transition, like in the standard mode-coupling theory for
systems with local dynamics only. The properties of the localized state are 
qualitatively similar to those of the localized state of a system with local 
dynamics only. The main difference is that at the same volume fraction the 
“cage size” is predicted to be larger in a system with particle size swaps 
than in the corresponding system without swaps. 

The properties of the localized state predicted by the present theory are 
qualitatively similar to those predicted by the replica approach <cit.>.
The main difference is that in our approach different “cage sizes” are 
obtained from the (approximate) calculation whereas the replica
ansatz assumes that there is a single “cage size” parameter. 

While solving numerically the time-dependent equations of motion for either
collective or tagged particle density correlations seems difficult, it might be 
possible to re-do the standard mode-coupling asymptotic analysis and
to check whether the power law approach to and departure from intermediate
time plateaus are still present, and to check whether there is an analogue
of the separation parameter determining the corresponding power laws exponents.
This interesting problem is left for a future study.



§ ACKNOWLEDGMENTS


I thank Th. Voigtmann for a reference to explicit expressions for binary hard
sphere mixture PY structure factors and E. Flenner for
comments on the manuscript.
I gratefully acknowledge partial support of NSF Grants No. DMR-1608086
and No. CHE-1800282.



§ 

99
BerthierBiroliRMP L. Berthier and G. Biroli,
Rev. Mod. Phys. 83, 587 (2011).
BerthierPNAS L. Berthier, P. Charbonneau, D. Coslovich, A. Ninarello,
M. Ozawa, and S. Yaida,
Proc. Natl. Acad. Sci. USA 114, 11356 (2017).
Grigera T. Grigera and G. Parisi, Phys. Rev. E 63, 045102(R) (2001).
Flenner2006 See, e.g., E. Flenner and G. Szamel, 
Phys. Rev. E 73, 061505 (2006).
BerthierCNO L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa,
Phys. Rev. Lett. 116, 238002 (2016).
NinarelloPRX A. Ninarello, L. Berthier, and D. Coslovich,
Phys. Rev. X 7, 021039 (2017).
OzawaPNAS M. Ozawa, L. Berthier, G. Biroli, A. Rosso and G. Tarjus,
PNAS 115, 6656 (2018).
WangNatComm L. Wang, A. Ninarello, P. Guan, L. Berthier, G. Szamel and
E. Flenner,
Nature Communications 10, 26 (2019).
Berthier2d L. Berthier, P. Charbonneau, A. Ninarello, M. Ozawa
and S. Yaida, arXiv:1805.09035.
WyartCates M. Wyart and M.E. Cates,
Phys. Rev. Lett. 119, 195501 (2017).
BBBT L. Berthier, G. Biroli, J.-P. Bouchaud, and G. Tarjus,
arXiv:1805.12378.
IkedaZamponiIkeda H. Ikeda, F. Zamponi, and A. Ikeda,
J. Chem. Phys. 147, 234506 (2017).
IkedaZamponi H. Ikeda and F. Zamponi,
arXiv:1812.08780.
BritoLernerWyart C. Brito, E. Lerner, and M. Wyart,
Phys. Rev. X 8, 031050 (2018).
Kapteijns 
G. Kapteijns, W. Ji, C. Brito, M. Wyart, and E. Lerner,
Phys. Rev. E 99 012106 (2019).
MDswap
L. Berthier, E. Flenner, C.J. Fullerton, C. Scalliet, and M. Singh,
arXiv:1811.12837.
Szamelswap G. Szamel,
Phys. Rev. E 98, 050601(R) (2018).
Goetzebook
W. Götze, Complex dynamics of glass-forming
liquids: A mode-coupling theory (Oxford University Press, Oxford, 2008).
HansenMcDonald J.P. Hansen and I.R. McDonald,
Theory of Simple Liquids
(Elsevier, Amsterdam, 2006).
CHess B. Cichocki and W. Hess, Physica A 141, 475 (1987).
KawasakiK. Kawasaki, Physica A 215, 61 (1995).
SL G. Szamel and H. Löwen, Phys. Rev. A 44, 8215 (1991).
PY1 J.L. Lebowitz and J. S. Rowlinson, J. Chem. Phys. 41, 133
(1964).
PY2 R. J. Baxter, J. Chem. Phys. 52, 4559 (1970).
Voigtmannthesis We used explicit formulas from
Th. Voigtmann, Ph.D. thesis, TU München, 2002.
GoetzeVoigtmann W. Götze and Th. Voigtmann,
Phys. Rev. E 67, 021502 (2003).
details We discretized the self-consistent equations (<ref>) using a 
grid of 401 equally spaced wavevectors, with q_mind_S=0 and 
q_maxd_S=80, where d_S is the diameter of the smaller sphere. 
The resulting equations were then solved by iteration, with integrations performed 
using a gaussian quadrature.


