Approach to classical liquids
based on functional renormalization group
Takeru Yokota

RIKEN iTHEMS

Refs.
TY, Haruyama, Sugino, Functional-renormalization-group approach to classical liquids with shortrange repulsion: A scheme without repulsive reference system, Phys. Rev. E 104, 014124 (2021)
TY, Physics-informed neural network for solving functional renormalization group on lattice,arXiv:2312.16038

Second Workshop on Fundamentals in Density Functional Theory (DFT2024), RIKEN Kobe campus, 22 Feb, 2024

About me
Main research topic: Functional renormalization group (FRG)
Graduate student at Kyoto U in 2014-2019

Prof. YoshidaDr. Naito

I started the study about
the combination of FRG and DFT.
w/ Prof. Kunihiro, Prof. Yoshida,
Dr. Naito, Mr. Kasuya

Researcher in the Sugino group at the Institute for
Solid State Physics, U Tokyo in 2020-2021

Prof. Sugino

Dr. Haruyama

I worked on FRG applications to
w/ Prof. Sugino, Dr. Haruyama
classical DFT.
Current position: Researcher in RIKEN iTHEMSMy recent interest is in machine-learning approaches to FRG.

About me
Main research topic: Functional renormalization group (FRG)
Graduate student at Kyoto U in 2014-2019

Prof. YoshidaDr. Naito

I started the study about
the combination of FRG and DFT.
w/ Prof. Kunihiro, Prof. Yoshida,
Dr. Naito, Mr. Kasuya

Today’s main topic

Researcher in the Sugino group at the Institute for
Solid State Physics, U Tokyo in 2020-2021

Prof. Sugino

Dr. Haruyama

I worked on FRG applications to
w/ Prof. Sugino, Dr. Haruyama
classical DFT.
Current position: Researcher in RIKEN iTHEMSMy recent interest is in machine-learning approaches to FRG.
Related topic (if time allows…)

Contents
Functional renormalization group (FRG) for classical liquids• FRG for classical scalar fields• Path integral for classical DFT• Modified FRG for classical liquids
Conventional approach: functional Taylor expansion• Functional Taylor expansion & hierarchical eq.
• Approach for hard-core case• Demonstration: 1D liquid
New approach: Machine learning• Physics-informed neural network for solving FRG• Demonstration: 0D O(N) model

Contents
Functional renormalization group (FRG) for classical liquids• FRG for classical scalar fields• Path integral for classical DFT• Modified FRG for classical liquids
Conventional approach: functional Taylor expansion• Functional Taylor expansion & hierarchical eq.
• Approach for hard-core case
TY, Haruyama, Sugino, PRE (2021)
Demonstration:
1D
liquid•New approach: Machine learning• Physics-informed neural network for solving FRG• Demonstration: 0D O(N) model

TY, arXiv:2312.16038

Contents
Functional renormalization group (FRG) for classical liquids• FRG for classical scalar fields• Path integral for classical DFT• Modified FRG for classical liquids
Conventional approach: functional Taylor expansion• Functional Taylor expansion & hierarchical eq.
• Approach for hard-core case• Demonstration: 1D liquid
New approach: Machine learning• Physics-informed neural network for solving FRG• Demonstration: 0D O(N) modelFunctional renormalization group (FRG) overview
Wegner, Houghton, PRA (1973), Wilson, Kogut,, PR (1974)
Polchinski, NPB (1984), Wetterich, PLB (1993)

FRG is…

Technique to analyze many-body systems• Versatile tool: condensed matter, statistical phys., high-energy phys.…
• Non-perturbative approach• Phase transition, phase diagram, thermodynamic quantities, correlations, …

Rigorous formulation of RG flow

• Derivation: path-integral formulation• Functional diﬀerential equation
Wetterich eq.

micro R

Wetterich, PLB (1993)∂k Rk
1∂kΓk[ϕ] = Tr
2 [ Γ(2)
]
[ϕ]
+
R
k
k
Governing eq. for “free energy” (I will explain it later)

Let’s see the formulation in the classical case.

gr

G

coa

ain rseingmacroClassical scalar field theory & effective action
Let us consider scalar fields ϕ(x) statistically fluctuating
following an action S[ϕ] (= Hamiltonian of the field βH[ϕ]).

1g22ϕ(x)( ∇ + m )ϕ(x) + ϕ(x)4E.g.) 1. Classical ϕ theory: S[ϕ] =)∫x ( 24!
2. Action for classical DFT ϕ(x) = ρ(x) (discussed later)
4

The partition function describes
the statistical properties:

Z[J] = 𝒟ϕe −S[ϕ]+ ∫x J(x)ϕ(x)∫Path integral= Summation over the configuration space(recall∑s1,…,sNe −βH({si}) in a spin system)

Classical scalar field theory & effective action
Let us consider scalar fields ϕ(x) statistically fluctuating
following an action S[ϕ] (= Hamiltonian of the field βH[ϕ]).

1g22ϕ(x)( ∇ + m )ϕ(x) + ϕ(x)4E.g.) 1. Classical ϕ theory: S[ϕ] =)∫x ( 24!
2. Action for classical DFT ϕ(x) = ρ(x) (discussed later)
4

The partition function describes
the statistical properties:

Z[J] = 𝒟ϕe −S[ϕ]+ ∫x J(x)ϕ(x)∫Similarly, the eﬀective action
Γ[ϕ] describes them.

Γ[ϕ] = supJ( ∫x

J(x)ϕ(x) − ln Z[J])• A kind of “free energy”• It provides a variational principle for ϕ (⟨ϕ⟩ satisfies δΓ[ϕ]/δϕ(x) = 0)
• Its derivatives give the correlation functions.

Classical scalar field theory & effective action
Let us consider scalar fields ϕ(x) statistically fluctuating
following an action S[ϕ] (= Hamiltonian of the field βH[ϕ]).

1g22ϕ(x)( ∇ + m )ϕ(x) + ϕ(x)4E.g.) 1. Classical ϕ theory: S[ϕ] =)∫x ( 24!
2. Action for classical DFT ϕ(x) = ρ(x) (discussed later)
4

Thedescribes
Wetterich’s formalism of−S[ϕ]+FRG ∫isbased on Γ[ϕ]The partition function
J(x)ϕ(x)xZ[J] = 𝒟ϕe
the statistical properties:∫
Wetterich, PLB (1993)

Similarly, the eﬀective action
Γ[ϕ] describes them.

Γ[ϕ] = supJ( ∫x

J(x)ϕ(x) − ln Z[J])• A kind of “free energy”• It provides a variational principle for ϕ (⟨ϕ⟩ satisfies δΓ[ϕ]/δϕ(x) = 0)
• Its derivatives give the correlation functions.

Wetterich’s formalism of FRG
In RG, the fluctuations of ϕ(x) are gradually taken fromhigh- to low-momentum components.
To implement this, we introduce a regulator term.
R̃k(p) Suppression ofSk[ϕ] = S[ϕ] + Sreg,k[ϕ]infrared modes

1
Sreg,k[ϕ] =ϕ(x)Rk(x − x′)ϕ(x′)
2 ∫x,x′
RG scale k

The regulator Rk is subject to1.
2.

lim R̃kUV(p) = ∞kUV→∞lim R̃kIR(p) = 0

kIR→0

Inclusion ofultraviolet
modes

p

Flow equation
Scale-dependent eﬀective action

Γ̃k[ϕ] = sup
J

( ∫x

)

J(x)ϕ(x) − ln Zk[J]Zk[J] = 𝒟ϕe −Sk[ϕ]+ ∫x J(x)ϕ(x)∫This interpolates S[ϕ] and Γ[ϕ]:For k → ∞, Γ̃k[ϕ] = S[ϕ] + Sreg,k[ϕ]∵ lim R̃k(p) = ∞For k → 0, Γ̃k[ϕ] = Γ[ϕ]∵ lim R̃k(p) = 0k→∞(saddle-point approx. in path int.)

k→0

The flow equation for Γ̃k[ϕ] is derived as follows…

Derivation of flow eq.
1∂k Γ̃k[ϕ] = −∂k Zk[Jsup,k]Zk[Jsup,k]

1=
𝒟ϕ∂k Sreg,k[ϕ]e −Sk[ϕ]+ ∫x Jsup,k(x)ϕ(x)Zk[Jsup,k] ∫1=∂k Rk(x − x′)
2 ∫x,x′∫ 𝒟ϕ′ϕ′(x)ϕ′(x′)e −Sk[ϕ′]+ ∫x Jsup,k(x)ϕ′(x)Zk[Jsup,k]δ 2 ln Zk[Jsup,k]1=∂k Rk(x − x′)+ ∂k Sreg,k[ϕ]2 ∫x,x′δJ(x)δJ(x′)

2pt propagatorδ ln Zk[Jsup,k]δ Γ̃k[ϕ]= ϕ(x)
Because= Jsup,k([ϕ], x) &
δJ(x)
δϕ(x)
2

We haveδ ln Z[Jsup,k]2−1

δ Γ̃k[ϕ]=(x, x′)
δJ(x)δJ(x′)( δϕδϕ )∫x′Functional inverse

A(x, x′)A −1(x′, x′′) = δ(x − x′′)

Exact equation for RG flow−1δ Γ̃k[ϕ]1∂kΓ̃k[ϕ] =∂k Rk(x − x′)(x′, x) + ∂kSreg,k[ϕ]2 ∫x,x′( δϕδϕ )Initial condition:

Γ̃kUV[ϕ] = S[ϕ] + Sreg,kUV[ϕ]Solution:Γ̃0[ϕ] = Γ[ϕ]Notation We use a diﬀerent notation from the Wetterich’s one,
where the flow equation is written as∂k Rk
1∂k Γk[ϕ] = Tr
2 [ Γ(2)]k [ϕ] + Rk

with the eﬀective average action Γk[ϕ] = Γ̃k[ϕ] − Sreg,k[ϕ].

Let us apply this formalism to classical liquids!

Classical liquids
Eve
ryw
her
aro
e
und
us!

Hamiltonian

HN =

N

2
pi

∑ 2m

+ VN({xi}i)i=1

N-body potential

VN ({xi}i) =∑v(xi, xj) +

i<j

∑

i<j<k

v3(xi, xj, xk) + ⋯Simple liquidsv(r)We only consider two-body interactions:

VN({xi}i) =∑i<j

Short-range
repulsionrLong-range
attractionv( | xi − xj | )Density field
We consider grand canonical ensemble (GCE).
Density (=positions of particles) is the fundamental degreeof freedom.

ϕ(x) = ρ(x) =

N̂

⟨∑
i=1

⟩

δ(x − x̂i)N,̂ x̂i are averaged following GCE.
GCE

*In quantum cases, Hohenberg-Kohn thm. is required for this statement.

“Density functional theory” naturally appears.
Let us find the action of density.

Partition function
Lue & Prausnitz, JCP (1998), Caillol, Mol. Phys. (2011)

To find the action, we consider the path integral of the
grand partition function∞

1−βHN({xi},{pi})+βμN
Ξ=e∑ N! ∫∫
p1,…,pN x1,…,xN
N=0
=

∞

N
n0

∑ N! ∫N=0−β( 12 ∫x,x′ ρN(x)v(x − x′)(ρN(x′) − δ(x − x′)))e
x1,…,xN
Density to be averaged ρN(x) =N∑δ(x − xi)i=1

n0 = Λ−3e βμ with de Broglie thermal wavelength Λ =2πβℏ2 /mΞ

Lue & Prausnitz, JCP (1998)∞n0N−β( 12 ∫x,x′ ρN (x)v(x − x′)(ρN (x′) − δ(x − x′)))=e∑ N! ∫
x1,…,xN
N=0=∫−β( 12 ∫x,x′ ρ(x)v(x − x′)(ρ(x′) − δ(x − x′)))
𝒟ρe

∞

n0N
×
δ(ρ(y) − ρN(y))∑ N! ∫∏
x1,…,xN y
N=0
= 𝒟ρe
∫

−β( 12 ∫x,x′ ρ(x)v(x − x′)(ρ(x′) − δ(x − x′)))

∞

n0N
×
𝒟ηe iβ ∫ dxη(x)(ρ(x)−ρN(x))∑ N! ∫∫
x1,…,xN
N=0=∫−β( 12 ∫x,x′ ρ(x)v(x − x′)(ρ(x′) − δ(x − x′)))
𝒟ρe

× 𝒟ηe
∫

∫x (iβη(x)ρ(x) + n0e −iβη(x))

Path integral
∫

𝒟ρ

∏

δ(ρ(y) − ρN(y)) = 1

yIntegral exp. of δδ(ρ( y) − ρN ( y)) = 𝒟ηe iβ ∫ dxη(x)(ρ(x)−ρN (x))∏∫
ySummation w.r.t. N
Particles
in GCE

density field ρmap+fictitious
potential ηAction for classical liquids

Ξ = 𝒟ρ e∫−S[ρ]The actionβS[ρ] = Sfree[ρ] +ρ(x)v(x − x′)(ρ(x′) − δ(x − x′))
2 ∫x,x′
Sfree[ρ] = ln

(∫

𝒟ηe

)

∫x (iβη(x)ρ(x) + n0e −iβη(x))

Action for classical liquids

Ξ = 𝒟ρ e∫−S[ρ]The actionβS[ρ] = Sfree[ρ] +ρ(x)v(x − x′)(ρ(x′) − δ(x − x′))
2 ∫x,x′
Sfree[ρ] = ln

(∫

𝒟ηe

)

∫x (iβη(x)ρ(x) + n0e −iβη(x))

Complicated!! !
Practically usable? I will come back to this point later.

Effective action
Ξ[U] = 𝒟ρe −S[ρ]+ ∫x βU(x)ρ(x)∫Grand partition function
w/ one-body ext. field U

Ω[U] = − β −1 ln Ξ[U]Grand potential

Effective action (=Free-energy density functional × β)Γ[ρ]( = βF[ρ]) = − supU[ ∫x

ρ(x)β(U(x) − μ) − βΩ[U]This comprehensively describes the system.

Variational principle: δF[ρ]/δρ(x) = μ
Thermodynamic quantities, correlations (derivatives)]

Problem! Regulator & initial condition
1
Sreg,k[ρ] =ρ(x)Rk(x − x′)ρ(x′)
2 ∫x,x′By choosing a regulator satisfying RkUV → ∞, we have theinitial condition ΓkUV[ρ] = S[ρ] + Sreg,kUV[ρ].

But the form of S[ρ] is complicated and does not look useful.
S[ρ] = ln

(∫

𝒟ηe

β
+
ρ(x)v(x − x′)(ρ(x′) − δ(x − x′))) 2 ∫x,x′∫x (iβη(x)ρ(x) + n0e −iβη(x))

How about discarding the condition RkUV → ∞ andintroducing another type of flow?

Deformation of the two-body potential
quadratic

1
1
S[ρ] + Sreg,k[ρ] = Sfree[ρ] +ρ(x)βv(x − x′)(ρ(x′) − δ(x − x′)) +ρ(x)Rk(x − x′)ρ(x′)
2 ∫x,x′
2 ∫x,x′Both the interaction and regulator terms are quadratic. Therefore, theregulator term can be regarded as the deformation of βv(x).
For more simplification,
we modify Sreg,k[ρ] as1
Sreg,k[ρ] =ρ(x)Rk(x − x′)(ρ(x′) − δ(x − x′))
2 ∫x,x′Then

1S[ρ] + Sreg,k[ρ] = Sfree[ρ] +ρ(x)βvk(x − x′)(ρ(x′) − δ(x − x′))
2 ∫x,x′with βvk(x) = Rk(x) + βv(x)Deformation of the two-body potential
We no longer consider the original RG flow. Thus we replace k with adeformation parameter λ ∈ [0,1].

1
S[ρ] + Sreg,λ[ρ] = Sfree[ρ] +ρ(x)βvλ(x − x′)(ρ(x′) − δ(x − x′))
2 ∫x,x′With the condition βvλ=1(x) = βv(x), the flow corresponds to thedeformation from a reference system with βv0 to the system of interest.

λ=0
βv0(x) = βvref(x)βvλ(x)λ=1
βv1(x) = βv(x)Flow eq.
The flow eq. is obtained in thesame manner as the Wetterich eq.

2−1δ Γ̃λ[ρ]1∂λΓ̃λ[ρ] =∂ R (x − x′)(x, x′) + ∂λSreg,λ[ρ]2 ∫x,x′ λ λ
δρδρ
(
)

Γ̃λ[ρ] = βFλ[ρ] (the free-energy density functional with vλ). Thus,2−1

δ βFλ[ρ]1∂λ βFλ[ρ] =∂λ βvλ(x − x′) ρ(x)ρ(x′) +(x, x′) − ρ(x)δ(x − x′)
2 ∫x,x′( δρδρ )Similar methods based on flow eqs.
Parola, Reatto, PRA (1985)
Hierarchical reference theory (HRT)• Almost the same equation is derived.
• vλ is decomposed as vλ = vrepulsive + vatt,λ• RG-like treatment of vatt,λ (inclusion from high- to low-momentum components ṽatt(p))• Phase transition (phase coexistence curve, critical exponent…)
Iso, Kawana, PTEP (2019)
Density RG• Flow eq associated with the spacial scale transformation
Polonyi, Sailer (2002), Schwenk, Polonyi (2004)
Quantum version: FRG-DFT• 1D nucleons [TY, Yoshida, Kunihiro PRC (2019)], 3D electrons [TY, Naito PRR (2021)]

How to solve functional differential eq.?
2−1

δ βFλ[ρ]1∂λ βFλ[ρ] =∂λ βvλ(x − x′) ρ(x)ρ(x′) +(x, x′) − ρ(x)δ(x − x′)
2 ∫x,x′( δρδρ )Initial condition: Fλ=0[ρ] = Fref[ρ]Solving functional diﬀerential eqs. is formidable!
Solution method?

• Conventional approximation: Taylor expansion(Next section)

TY, Haruyama, Sugino, PRE (2021)• New approach: Machine learning(Next to next section)

TY, arXiv:2312.16038

Contents
Functional renormalization group (FRG) for classical liquids• FRG for classical scalar fields• Path integral for classical DFT• Modified FRG for classical liquids
Conventional approach: functional Taylor expansion• Functional Taylor expansion & hierarchical eq.
• Approach for hard-core case
TY, Haruyama, Sugino, PRE (2021)
Demonstration:
1D
liquid•New approach: Machine learning• Physics-informed neural network for solving FRG• Demonstration: 0D O(N) modelFunctional Taylor expansion
Liquids are almost homogeneous. We expand Fλ[ρ] around ρhom.
2−1

δ βFλ[ρ]1∂λ βFλ[ρ] =∂λ βvλ(x − x′) ρ(x)ρ(x′) +(x, x′) − ρ(x)δ(x − x′)
2 ∫x,x′( δρδρ )(1)βFλ[ρ] = βFλ[ρhom] + dxβFλ ([ρhom], x)(ρ(x) − ρhom) + ⋯∫Hierarchical eq. for Fλ[ρhom], F (1)[ρhom], F (2)[ρhom], …
λλThis should be truncated at some order tofacilitate numerical computations.

Hierarchical eq. for correlation functions
Fλ[ρhom], Fλ(1)[ρhom], Fλ(2)[ρhom], … are related to the free energy,

chemical potential, and distribution functions through variational eq.
& Ornstein-Zernike eq.
Hierarchical eq. for distribution functions
βFλ
ρhom
∂λ
=−
∂λ βvλ(r)gλ(2)(r)N
2 ∫r

0th order
Free energy

1st order2∂λ βμλ =Chemical potential

ρ
(3)
(2)
∫
ρhom ∫r ∂λ βvλ(r)gλ(2)(r) + hom
∂
βv
(r
−
r′
)
g
(r,
r′
)
−
g
(r − r′)]
[
λ
λ
λ
λ
2
r,r′

1 + ρhom ∫r [gλ(2)(r) − 1]nn( ≥ 2)th order ∂λgλ (r1, ⋯, rn−1) = − ∑ ∂λ [βvλ(ri − rj)] gλ (r1, ⋯, rn−1)(n)n-particle

distribution function

(n)

i<j

+ngλ(n)(r1, ⋯, rn−1)∂λ μλ + ρhom∂λ μλ
−ρhom

n

∫r ∑
i=1

∫r

(n+1)
(n)
g
(r
,
⋯,
r
,
r)
−
g
(r1, ⋯, rn−1)]
1
n
[ λ
λ

∂λ [βvλ(r − ri)] gλ(n+1)(r, r1, ⋯, rn−1)

2
ρhom
−
∂λ [βvλ(r − r′)] [gλ(n+2)(r, r′, r1, ⋯, rn−1) − gλ(2)(r − r′)gλ(n)(r1, ⋯, rn−1)]2 ∫r,r′Hard-core divergence?
v(r)

βFλ
ρhom
∂λ
=−
∂λ βvλ(r)gλ(2)(r)N
2 ∫r

Short-range
repulsionr

At first glance, the repulsive coreleads to a divergence.

Long-range
attraction

In HRT, the repulsive core is treated as the reference in order
to avoid to treat it in the flow eq.

vλ = vref + vatt,λ∂λvλ = ∂λvatt,λ(2)However, this divergence does not occur since g (r) is
λsuppressed.

Cavity distribution function
In fact, the divergence is removed when we introduce

the cavity distribution function

Exact eqs. for cavity distribution functions0thnβ ∑i<j vλ(ri−rj) (n)(n)
yλ (r, ⋯, rn−1) = egλ (r, ⋯, rn−1)

βFλ ρhom∂λ=∂λ fλ(r)yλ(2)(r)N
2 ∫r

No divergence!

(Only Mayer function fλ(r, r′) = e −βvλ(r,r′) − 1 appears)2ρ(2)−βvλ(r)−βvλ(r′) (3)∫−ρhom ∫r ∂λ fλ(r)yλ(2)(r) − hom
∂
f
(r
−
r′
)
e
y
(r,
r′
)
−
y
[
λ
λ
λ
λ (r − r′)]
2
r,r′

1st

∂λ βμλ =nth∂λ yλ(n)(r1, ⋯, rn−1) = ∂λ βμλ nyλ(n)(r1, ⋯, rn−1) + ρhom [e −β ∑i∫r(

1 + ρhom ∫r [e −βvλ(r)yλ(2)(r) − 1]
n−1

+ρhom

∫r

∂λ [e

n−1

−β ∑i

] yλ

vλ(ri−r)−βvλ(r)Method w/o
hard-core reference!

vλ(ri−r)−βvλ(r) (n+1)
yλ (r1, ⋯, rn−1, r) − yλ(n)(r1, ⋯, rn−1)]

)

(n+1)

(r, r1, ⋯, rn−1)

2
n−1
n−1
ρhom
+
∂λ fλ(r − r′)[e −β ∑i vλ(ri−r)−β ∑i vλ(ri−r′)−βvλ(r)−βvλ(r′)yλ(n+2)(r, r′, r1, ⋯, rn−1) − yλ(2)(r − r′)yλ(n)(r1, ⋯, rn−1)]2 ∫r,r′

The repulsive reference is unnecessary.
We use the free-gas reference v0 = 0.

βF0[ρ] =∫x

ρ(x)(ln (Λ3ρ(x)) − 1)

Truncation up to 2nd order
0th orderβFλ ρhom∂λ=∂λ fλ(r)yλ(2)(r)N
2 ∫r
1st order∂λ βμλ =−ρhom ∫r ∂λ fλ(r)yλ(2)(r) −

2
ρhom

2

∫r,r′ ∂λ fλ(r − r′)[e −βvλ(r)−βvλ(r′)yλ(3)(r, r′) − yλ(2)(r − r′)]1 + ρhom ∫r [e −βvλ(r)yλ(2)(r) − 1]

2nd order∂λ yλ(2)(r1) = ∂λ βμλ 2yλ(2)(r1) + ρhom [e −βvλ(r1−r)−βvλ(r)yλ(3)(r1, r) − yλ(2)(r1)]
∫r
(
)
+ρhom

∫r

∂λ [e

−βvλ(r1−r)−βvλ(r)(3)] yλ (r, r1)

Approximation is needed

2
ρhom
+
∂λ fλ(r − r′)[e −βvλ(r1−r)−βvλ(r1−r′)−βvλ(r)−βvλ(r′)yλ(4)(r, r′, r1) − yλ(2)(r − r′)yλ(2)(r1)]2 ∫r,r′Truncation up to 2nd order
0th order

We close the hierarchy using
Kirkwood superposition approximation (KSA)

Kirkwood, JCP (1935)• Ignorance of mean many-body potential
βFλ ρhom(2)
Exact for low densities
•
∂λ
=
∂
f
(r)y
(r)
(3)
(2)
(2)
(2)
λ λr ) ≈1 λ yλ (r1)yλ (r)yλ (r − r1)N
2 ∫yλr (r,1st order∂λ βμλ =yλ(4)(r, r′, r1) ≈ yλ(2)(r)yλ(2)(r′)yλ(2)(r1)yλ(2)(r − r′)yλ(2)(r′ − r1)yλ(2)(r1 − r)−ρhom ∫r ∂λ fλ(r)yλ(2)(r) −

2
ρhom

2

∫r,r′ ∂λ fλ(r − r′)[e −βvλ(r)−βvλ(r′)yλ(3)(r, r′) − yλ(2)(r − r′)]1 + ρhom ∫r [e −βvλ(r)yλ(2)(r) − 1]

2nd order∂λ yλ(2)(r1) = ∂λ βμλ 2yλ(2)(r1) + ρhom [e −βvλ(r1−r)−βvλ(r)yλ(3)(r1, r) − yλ(2)(r1)]
∫r
(
)
+ρhom

∫r

∂λ [e

−βvλ(r1−r)−βvλ(r)(3)] yλ (r, r1)

Approximation is needed

2
ρhom
+
∂λ fλ(r − r′)[e −βvλ(r1−r)−βvλ(r1−r′)−βvλ(r)−βvλ(r′)yλ(4)(r, r′, r1) − yλ(2)(r − r′)yλ(2)(r1)]2 ∫r,r′1D solvable hard rod
Let us apply our method to thederivation of the pair distributionfunction (PDF) & free energy in a
1D solvable model.

Model

Exclusion by surrounding particles(shell structure)…
Repulsivecore
Yarnell, Katz, Wenzel, Koenig PRA (1973)

Archer, Chacko, Evans, J. Chem. Phys. (2017)βv(x)βv(x) =σ−zpσppair distribution function

Typical behavior of PDF (Ar)σp + σx∞( | x | ≤ σ)
−zp(σ + σp − | x | ) (σ < | x | ≤ σ + σp)0(σ + σp ≤ | x | )In the case of σp < σ, thermodynamic
quantities & pair correlation function are
obtained exactly.
Takahashi, Proc. Phys.-Math. Soc. Jpn. (1942)
Percus, J. Stat. Phys. (1982)
Brader, Evans, Physica A (2002)

Choice of evolution

Naive adiabatic connection cannot be used for hard corefλ(x) = e −λβv(r) − 1 suddenly changes at λ = 0Instead, we introduce the evolution of hard-core diameter∞hardβvλ (x) ={0( | x | ≤ λσ)(λσ < | x | )β v(x)

λσ

xThis is compatible with KSA.

ρhomσ is the only dimensionless parameter

Small σ corresponds to small ρhom

At the beginning of the flow (small σ), KSA is very accurate

The evolution of the attractive part is the adiabatic connection.

Benchmarks
We calculate the free energy & pair distribution function.

• Runge-Kutta & numerical spatial integral
Benchmarks• Exact result• Integral-equation method

A widely used method for classical liquids.

Ornstein-Zernike eq. (Schwinger-Dyson eq.)g(x) = 1 + c(x) + ρ0 dx′c( | x − x′| )[g(x′) − 1]∫c(x − x′) = δ(x − x′) −+Closure (approximation)g(x) ≈ g(x, c(x))δβF[ρ]
δρ(x)δρ(x′)

Closureg(x) = exp [−βv(x) + g(x) − c(x) − 1]Hypernetted chain (HNC)• Based on 2nd order functional Taylor expansion for βF[ρ]& g (2)(x, x′) = ρ0−1ρ(x | x′) (Percus’ test particle method)g(x) = exp [−βv(x)] [g(x) − c(x)]Percus-Yevick equation (PY)• Linearization of HNC• Exact for one-dimensional hard rod w/o attractive forces

Free energy for hard rod (zpσ = 0)β v(x)

Excess free energy
σ(hypernetted chain)(=PY)

Vertex expansion(VE) shows moreaccurate result than
HNC

ρhomσ

x

Distribution functions for Hard rod (zpσ = 0)
Two-particle distribution function ρhomσ = 0.55

FRG-DFTFRG

FRG-DFT
FRG
HNC
Exact (=PY)
(PY)

HNC
Exact (PY)
(=PY)

1.5
g(x)

2

c.f.) lower density ρhomσ = 0.4

1.0

g(x)

0.5
0.0

0

1

2

3

4

5

x/æ

1

0

Height of 1st peak&
position of 2nd peak

0

1

2

3

4

5

are described by VE moreaccurately than HNC

x/æ
KSA + diameter evolution successfully describe hard part!

Free energy with attractive force (zpσ = 1, σp /σ = 0.9)
β v(x)

Excess free energy

0.5ØF ex/N

0.4
0.3

σ

exact
FRG
flow
HNC (hypernetted chain)
PY (Percus-Yevick eq.)

−zpσp

PY deviates from exact
result due to theattractive force

0.2

VE shows moreaccurate result than
HNC and PY

0.1
0.0
0.0

0.2
ρhomΩæσ0.4

σp + σ xDistribution functions (zpσ = 1, σp /σ = 0.9)
Two-particle distribution function

ρhomσ = 0.55
FRG
Flow

g(x)

3

HNC
PY
Exact

c.f.) lower density ρhomσ = 0.4
FRG
Flow
FRG
HNC
PY
Exact

g(x)

2

1

2
0

0

1

2

3

4

5

x/æ

1
0

Height of 1st peak
Is described by VE moreaccurately than HNC and PY

0

1

2

Diﬀerence from hard rod:
Height of 1st peak increases
due to attractive force

3

4

5

x/æ

VE accurately incorporates not only the hard part but also the attractive part!

Contents
Functional renormalization group (FRG) for classical liquids• FRG for classical scalar fields• Path integral for classical DFT• Modified FRG for classical liquids
Conventional approach: functional Taylor expansion• Functional Taylor expansion & hierarchical eq.
• Approach for hard-core case• Demonstration: 1D liquid
New approach: Machine learning• Physics-informed neural network for solving FRG• Demonstration: 0D O(N) model

TY, arXiv:2312.16038

Taylor-expansion-based method and its limitations
Let us switch the notation to Wetterich’s original one.
2

−1

δ Γk
1
∂k Γk[φ] = Tr ∂k Rk[φ] + Rk2
( δφδφ)Γk[φ] = Γ̃k[φ] − Sreg,k[φ]Tr: trace w.r.t. coordinates& internal DOF

We have applied the functional Taylor expansion to the flow equation.
But, there are limitations…
❌ The solution is applicable to limited configurations of φ(x)• φ(x) ≈ φtarget(x)• High cost for the analysis of complicated field configurations, such as crystal.
❌ Usually, improving the order of the truncation is not easy.

• Many coordinate (momentum) integrals

Other solution methods?

Another attempt: FDE as high-dim. PDE
FDE: Infinite-dimensional partial diﬀerential equation (PDE)
To realize numerical analysis, the input dimensions should be truncated.
Method 1) Introduction of finite spatial lattice
φ = {φn,α}n,α, n = (n1, …, nd), 0 ≤ ni < L, α = 1,…, NIDOF (internal DOF)
Total DOF of φ: NDOF = L d NIDOF

Method 2) Basis function expansion
φ(x, α) ≈NDOF∑

See, e.g., Venturi, PR (2018),
Venturi, Dektor, Res. Math. Sci (2021)

φibi(x, α) = φ ⋅ b(x, α) (bi(x, α): orthonormal basis function)

i=1

2

−1

δ Γk
1
∂k Γk[φ] = Tr ∂k Rk[φ] + Rk2
( δφδφ)(NDOF + 1)-dim. PDE

2

−1

∂ Γk
1
∂k Γk(φ) = tr ∂k Rk(φ) + Rk2
( ∂φ∂φ)The original continuum theory is obtained for increasing NDOF.

Failure of computational grids
How can we solve (NDOF + 1)-dim. PDE with large NDOF

Taking computational grids for φ is infeasible.
# of grid points ∼ exp(NDOF)

if grids are assigned for each φn,αGrid-based methods (finite-element method, Runge-Kutta, …) can not be used.

c.f.) Some people attempt to use tensor decomposition to mitigate
computational complexity of representing the solution for linear
FDE, such as Hopf eq.
But only the results with NDOF ≲ 6 has been reported…
Venturi, PR (2018), Venturi, Dektor, Res. Math. Sci (2021)

Grid-free method for high-dim. PDE is required

Physics-informed neural network (PINN)
ML-based grid-free approach: PINN

…
…

…

…

…

activationx ∈ Ω (Domain)x ∈ ∂Ω

…

xN…
…

f(x) ≈ fNN(x; θ)

x1
x2
…

NN assumption
for solution

IDE( f(x), ∂μ f(x), ⋯) = 0
IBC( f(x), ∂μ f(x), ⋯) = 0
activation

PDE:
Boundary condition:

I.E. Lagaris, A. Likas, D. I. Fotiadis, IEEE Transactions on Neural Networks (1998)
M. Raissi, P. Perdikaris, G. E. Karniadakis (2017)
M. Raissi, P. Perdikaris, G.E. Karniadakis, Journal of Computational Physics (2019)fNN(x; θ)

Learning parameter θ

Loss function (taking the minimum when PDE & BC are satisfied)
例: L(θ) = average∥IDE( fNN(x; θ), ∂μ fNN(x; θ), ⋯)∥2 + average∥IBC( fNN(x; θ), ∂μ fNN(x; θ), ⋯)∥2
x∈Ω

x∈∂Ω

PINN’s capability of handling high-dim. PDEs
Advantage of PINN: Grid free

Applicable to high-dim. PDEs

Applications to high-dim. PDEs
• L. Guo, H. Wu, X. Yu, T. Zhou, Computer Methods in Applied Mechanics and Engineering (2022)• D. He, S. Li, W. Shi, X. Gao, J. Zhang, J. Bian, L. Wang, T.-Y. Liu, In International Conference on Artificial
Intelligence and Statistics (2023)• J. Cen, X. Chen, M. Xu, Q. Zou, arXiv:2305.06863.
• K. Tang, X. Wan, C. Yang, Journal of Computational Physics, 476 (2023)• Z. Hu, K.Shukla, G. E. Karniadakis, K. Kawaguchi, arXiv:2307.12306• Z. Hu, Z. Yang, Y. Wang, G. E. Karniadakis, K. Kawaguchi, arXiv:2311.15283• Z. Hu, Z. Shi, G. E. Karniadakis, K. Kawaguchi, arXiv:2312.14499

105-dim. PDE

Hamilton-Jacobi-Bellman eq.
Black-Scholes-Barenblatt eq.

Moreover, the solution is obtained
simultaneously for a domain of the
inputs rather than one input.

PINN for Wetterich equation
PINN allows us to solve the Wetterich equation with large NDOF!

(k, φ)

NN(θ)Γ θk (φ)

Lθ = LθPDE + LθBC
LθPDE = 𝔼φ∼𝒫φk ∼ 𝒫k

LθBC =

1θ∂k Γk (φ) − tr
2

∂k Rk

∂2Γθk (φ)+ Rk( ∂φ∂φ
)

−1

2

θ
Γ
(φ) − S(φ))
k(UV]φ ∼ 𝒫φ [𝔼

2• The expectation is evaluated on a finite number of collocation points (k, φ).
• 𝒫φ,k is some probability distribution.

PINN for Wetterich equation
PINN allows us to solve the Wetterich equation with large NDOF!

(k, φ)

NN(θ)Γ θk (φ)

Lθ = LθPDE + LθBC
LθPDE = 𝔼φ∼𝒫φk ∼ 𝒫k

LθBC =The initial condition may be implemented
directly in NN.
Hereafter, we omit LθBC.
Lagaris, Likas, Fotiadis, IEEE Transactions on Neural Networks (1998)1θ∂k Γk (φ) − tr
2

∂k Rk

∂2Γθk (φ)+ Rk( ∂φ∂φ
)

−1

2

θ
Γ
(φ) − S(φ))
k(UV]φ ∼ 𝒫φ [𝔼

2• The expectation is evaluated on a finite number of collocation points (k, φ).
• 𝒫φ,k is some probability distribution.

0-dim. O(N) model• NDOF = N

1 2 2 g 22
S(φ) = m φ + (φ )2
4!

The Wetterich equation is an (N + 1)-dim. PDE.

We can investigate the scalability with NDOF by increasing N.
* We do not reduce Wetterich eq. to 2-dim. PDE with k and ρ = φ 2 /2

• Exact results and results by perturbative, large-N expansions areavailable.

E.g., Keitel, Bartosch, JPA (2012)

The perturbative region is given by g̃ = Ng/m 4 ≪ 1

•
2
γ(k,
φ)
σ
=
∂
We
calculateandself-energy
φγ(k, φ).
•2e −2l, l = ln(kUV /k)Regulator: Rkαα′ = kUV

•
2 2
Mass
squared:m/kUV = 0.01 ( ≪ 1 to validate the UV saddle-point cond.)•

Neural network for effective action
We replace the “nontrivial part” of the eﬀective action with an NN.
Interaction-induced RG partθΓk (φ) = S(φ) + ΔSfree(k) + γθ(k, φ) (“nontrivial part”)
Constant shiftassociated with free part∂2Sfree(φ)
1∂kΔSfree(k) = tr ∂k Rk+ Rk
2( ∂φ∂φ)This satisfies γθ(kUV, φ) = 0−1

NN ansatz subject to γθ(kUV, φ) = 0 (hard implementation of the init. cond.)γθ(k, φ) ≈ NNθ(l, φ) − NNθ(0,φ)…

inputs…

Field φN+1
inputs…
…

NNθ

3 hidden layers
256 units/layer
Softplus activation⊖ γ(l, φ)

…
…

l=0

…
…

…

• 3 hidden layers• 256 units/layer• Diﬀerentiable softplus activation

RG scale l

N+1…

NNθ(l, φ)l = ln(kUV /k)NNθPretrainingLθ =𝔼

φ ∼ 𝒫φ1θ∂lΓl (φ) − tr2l ∼ 𝒫l∂l Rl∂2Γθl (φ)( ∂φ∂φ

−1

+ Rl

2

)

The matrix must be regular during the training.
In our experience, this is frequently broken with randomly chosen θ.

Pretraining with some approximate analytic
results remedies this problem
We use 1st-order perturbative result:(γ(l, φ; θ) − γ (l, φ)) ][φpreLθ = 𝔼φ∼𝒫l ∼ 𝒫l1pt

2Other details of numerical procedure
−1

Lθ =

𝔼

φ ∼ 𝒫φl ∼ 𝒫l22 θ∂
Γl (φ)
1θ∂lΓl (φ) − tr ∂l Rl
+ Rl
2
∂φ∂φ
(
)

• 𝒫l: uniform distribution in [0,lend] with lend = 5
2:𝒫issampledfollowing∥φ∥N(0,N/m) (w/o sign)• φn̂ = φ/∥φ∥ is uniformly sampled* Other choices such as N(0, m −21) fail to sample the neighborhoods of
φ = 0 due to curse of dimensionality• 500 collocation points are used to evaluate the expectation.
• Adam optimizer• Pytorch• The matrix inverse is evaluated by direct method• This may not be eﬃcient but is easy to implement (torch.linalg.inv)• More eﬃcient way: Hutchinson trace estimator (future work) Hutchinson, Simul. Comp. (1990)

Code: https://github.com/TakeruYokota/PINNLFRGComputational time & convergence
We conduct computations for all the combinations ofg̃ = Ng/m 4
N = 1,10,100 and g̃ = 0.1,1,10

Learning curve & histories of physical quantities(N = 100 & g̃ = 1 case)Comp. time on NVIDIA A100 GPU

0.15
0.5

°0.5

0.14

0

°1.0

3.6

°1.5

¢æ/m2

3.8

æ/m2

Lµ (logscale)

0

4.0

0.0

∞

•
•

Learning rate (Wetterich): 10−4 with exponential decay factor 0.99999
Learning rate (pretraining): 10−3

0.13

0.12

0

3.4

°2.0

0.11

°2.5

0

3.2

°3.0
0.0

0.5

Iteration

1.0
£106

0.0

0.5

Iteration

1.0
£106

0

0.0

0.5

Iteration

1.0
£106

0

0.00

æ(l, 0)/m2

N = 1 and g̃ = 1 case
0.2

0.0

RG scale l dependence at φ = 0

φ dependencel at l = lend

0

1

2

3

4

5

(b)

(a)

0.05

∞(lend , ')PINN-LFRG (our method)
Pretrained model
Exact
Perturbation
Large-N

0.10∞(l, 0)

0.4

0.3
0.2
0.1

æ(lend , ')/m2

æ(l, 0)/m2

0.00
0.4

0.2

0.0
0

1

2

3

l

∞(lend , ')

(b)

4

5

l = ln(kUV /k)

0.4

0.2
°1.00 °0.75 °0.50 °0.25

0.25

0.50

0.75

m'

• γ(l, φ) is simultaneously obtained for a domain of φ in our method

0.3(PINN-LFRG).

0.2• PINN-LFRG shows accurate results compared to 1st-order

0.1

perturbation & leading-order large-N expansion

')/m2

0.00

0.4

1.00

æÆ (l, 0)/m

N = 100 and g̃ = 1 case
0.10
0.05
0.00

RG scale l dependence at φ = 0(a)PINN-LFRG
Exact
Perturbation
Large-N2(b)∞(lend , 'eÆ )∞(l, 0)4

φ dependencel at l = lend
0

1

2

3

4

5

4.25
4.00

Æ = 1, . . . , 100 (PINN-LFRG)

3.75

æÆ (lend , 'eÆ )/m2

æÆ (l, 0)/m2

0Æ = 1, . . . , 100 (PINN-LFRG)

0.15
0.10
0.05
0.00
0

1

2

3

l

5

0.16
Æ = 1, . . . , 100 (PINN-LFRG)
0.15
0.14
°1.00 °0.75 °0.50 °0.25∂2Reduced self-energy σα(l, φ) =γ(l, φ)∂φα2(b)∞(lend , 'eÆ )

4

0.17

4.25

0.00

0.25

0.50

0.75

m'

• Except for γ(l, 0), PINN-LFRG results are given by N = 100 lines
4.00

Æ = 1, . . . , 100 (PINN-LFRG)

3.75

corresponding to the N-direction in φ space.

æÆ (lend , 'eÆ )/m2

• PINN-LFRG shows comparable results with large-N expansion,
0.17
0.16

which should beÆ accurate
for N = 100.
= 1, . . . , 100 (PINN-LFRG)• O(N) symmetry is reproduced in PINN-LFRG.
0.15
0.14

°1.00 °0.75 °0.50 °0.25

0.00

m'

0.25

0.50

0.75

1.00

1.00

Results at φ = 0 for different N and g̃
Relative error compared to exact results (minus indicates underestimation)
γ = γ(lend, 0)1 N
σ=σα(lend, 0)N∑
α=1Absolute relative error of ∞(lend , 0) (%)Δσ =PINN-LFRG (g̃ = 0.1)
PINN-LFRG (g̃ = 1)

300

PINN-LFRG (g̃ = 10)
Perturbation (g̃ = 0.1)

250

Perturbation (g̃ = 1)

1 N
2(σα (lend, 0) − σ)∑N α=1

• For all g̃ and N, the errors of PINN-LFRG are within 3% for γ and 1% for σ.

Perturbation (g̃ = 10)200

Large-N (g̃ = 0.1)

• Even when 1/N and g̃ are not small,

Large-N (g̃ = 1)
150

Large-N (g̃ = 10)

PINN-LFRG shows accurate results.

100
50
0
100

101

N

102

NNs are promising approximations
independent of the existence of asmall parameter.

Summary
I showeda functional renormalization group (FRG) formalism for classical liquids.

=

FRG for classical DFT• FRG for classical scalar fields• Action for classical DFT• Modified FRG for classical liquids
Two approaches for solving FRG

Functional Taylor expansion

TY, Haruyama, Sugino, PRE (2021)• Solution method with no need for a hard-core reference• Demonstration in 1D liquids.
Physics-informed neural network solver• FRG eq: a kind of high-dimensional PDE• PINN: a cutting-edge solver for high-D PDE• Demonstration in the 0D scalar model
Outlook• 3D liquids• PINN for classical DFT

TY, arXiv:2312.16038

Backup

Kirkwood superposition approximation (KSA)
Free particle’s distribution in a potential U(x)

ρ(x) = ρ0e −βU(x)xx′Pair distribution= particle distribution around a fixed particle

fixedW(x − x′)Mean force potential (MFP)
generated by the fixed& surrounding particles

g (2)(x, x′) = ρ0−1ρ(x | x′) = e −βW(x−x′)N-particle distribution= (particle distribution around N − 1 fixed particle) x (N − 1)-particle distribution

g (3)(x1, x2, x3) = ρ0−1ρ(x1 | x2, x3)g (2)(x2, x3) = e −βW(x1|x2,x3)e −βW(x2−x3)

KSA = Superposition of MFPs

W(x1 | x2, x3) ≈ W(x1 | x2) + W(x1 | x3)

Removal of mean many-body forcesg (3)(x1, x2, x3) ≈ g (2)(x1, x2)g (2)(x2, x3)g (2)(x3, x1)

4-body W(x1 | x2, x3, x4) ≈ W(x1 | x2) + W(x1 | x3) + W(x1 | x4)
g (4)(x1, x2, x3) = ρ0−1ρ(x1 | x2, x3, x4)g (3)(x2, x3, x4) = e −βW(x1|x2,x3,x4)g (3)(x2, x3, x4)= g (2)(x1, x2)g (2)(x1, x3)g (2)(x1, x4)g (2)(x2, x3)g (2)(x1, x4)g (2)(x3, x4)

Hypernetted chain
Percus’s test particle methodxn(x | x′) = ρ0g (2)(x, x′)i.e., g (2)( ⋅ , x′) is obtained as the density in thepresence of the external field Ux′(x) = v(x − x′)

x′
fixed

Ux′(x) = v(x − x′)

Variational eq.
δβF[ρ0g (2)( ⋅ , x′)] = β(μ − v(x − x′))δρ(x)Functional Taylor expansion of the excess part around ρ(x) = ρ0βF[ρ] ≈

∫x

ρ(x)(ln ρ(x)/ρ0 − 1) + βFex[ρ0] +∫x

βμ(ρ(x) − ρ0) −1
c(x − x′)(ρ(x) − ρ0)(ρ(x′) − ρ0)2 ∫x,x′βF[ρ] = βFfree[ρ] + βFex[ρ]ClosureFrom variational eq. & Taylor exp.,
g (2)(x) = exp −βv(x) + ρ0 c(x − x′′)(g (2)(x, x′′) − 1)
∫x′′
[
]

OZ eq.

g (2)(x) = exp [−βv(x) + g (2)(x) − c(x) − 1]Machine learning for partial differential equations (PDE)
Recently, there have been many applications of machine learning to PDEs

• Physics-informed neural network (PINN)
Psichogios, Ungar, AIChE (1992)
Lagaris, Likas, Fotiadis, IEEE Transactions on Neural Networks (1998)
Raissi, Perdikaris, Karniadakis (2017)
Raissi, Perdikaris, Karniadakis, Journal of Computational Physics (2019),…

• Rayleigh-Ritz variational method
E, Yu, Comm. Math. Stat. (2018)
Khoo, Lu, Ying, Res. Math. Sci. (2018),…

• Backward stochastic diﬀerential equation
E, Han, Jentzen, Comm. Math. Stat. (2017)
Han, Jetzen, E (2018)
Rassi (2018)
Beck, E, Jentzen, J. Nonlinear Science (2019),…

Applicable to various
types of PDEs…