Finding flux vacua

Finding flux vacua#

What’s in this notebook? This notebook demonstrates how to use numerical optimisers to construct flux vacua by either finding solutions to the \(F\)-term conditions \(D_IW=0\) or to the extremum conditions \(\partial_I V=0\) (with \(D_IW\neq 0\)).

In this notebook, you will learn:

How to set up a flux compactification model with FluxEFT

How to find SUSY flux vacua by solving \(D_aW = 0\) with newton_method_flux_vacua

How to verify solutions by checking the Hodge decomposition of \(G_3\) and the F-term residual

How to find non-SUSY vacua by minimising the scalar potential \(V\)

Prerequisites: NB02: JAXVacua overview.

(Created: Andreas Schachner, June 25, 2024)

Setup#

# General imports
import warnings
import time
import numpy as np
from tqdm.auto import tqdm
from scipy.optimize import root

# JAX imports
from jax import vmap
import jax
import jax.numpy as jnp
jax.config.update("jax_enable_x64", True)

# Plotting
import seaborn as sn
import matplotlib.pyplot as plt
cmap = sn.color_palette("viridis", as_cmap=True)

# JAXVacua
import jaxvacua as jvc

warnings.filterwarnings('ignore')

Finding flux vacua at \(h^{1,2}=2\)#

We look at the typical example of the degree 18 hypersurface in \(\mathbb{CP}[1,1,1,6,9]\)

h12 = 2
model = jvc.FluxEFT(h12=h12,model_ID=1)

Solutions to the \(F\)-term conditions#

First, let us look at the following solution

# Moduli values
z0 = jnp.array([0.31686516+3.2301486j , 0.37113099+3.07462404j])

# Axio-dilaton values
tau0 = -0.29834961+6.80357615j

# Choices of fluxes
fluxes = jnp.array([4.,  -3.,  -2.,  -2.,   3.,   2.,  39., -13.,  -4.,   0.,  -0.,  0.])

This is a solution to the F-term conditions \(D_IW=0\) as can be verified by computing

model.DW(z0,jnp.conj(z0),tau0,jnp.conj(tau0),fluxes)

Of course, the numerical tolerance \(|D_IW|<\epsilon\) is only as good as the precision used above.

We now select some starting point and try to find the above solution through numerical minimisation.

# Set starting point
moduli = 2*jnp.ones(h12)*1j
tau = 1j*2.

First, let us try Newton’s method

# Set parameters for minimsation
step_size_Newton = 1
tol = 1e-10

print("Jit compiling...")
# Compute objective
DW = model.DW(moduli,jnp.conj(moduli),tau,jnp.conj(tau),fluxes)
# Compute Jacobian
dDW = model.dDW_real(moduli,tau,fluxes)

print("Start minimisation...")
counter = 0
tic = time.time()
res = 10.
while res>tol:
    
    # Compute objective
    DW = model.DW(moduli,jnp.conj(moduli),tau,jnp.conj(tau),fluxes)
    
    # Reshape DW to be a real vector
    DW_real = jnp.concatenate([jnp.real(DW),jnp.imag(DW)])
    
    # Compute Jacobian
    dDW = model.dDW_real(moduli,tau,fluxes)
    
    # Compute Newton step
    delta = -jnp.linalg.inv(dDW)@DW_real
    
    # Taking step in moduli direction
    moduli = moduli + step_size_Newton*(delta[h12+1:2*h12+1]+1j*delta[0:h12])
    
    # Taking step in tau direction
    tau = tau + step_size_Newton*(delta[2*h12+1]+1j*delta[h12])
    
    # Compute new residual
    res = jnp.sum(jnp.abs(DW_real))
    counter+=1
    
    print(f"residual: {res}                     counter: {counter}           ", end="\r")
    
    if counter>10**3:
        print("")
        print(f"STOPPED RUN: residual: {res}                     counter: {counter}           ")
        break
        
toc = time.time()
print("")
print(f"Finished in {np.around((toc-tic)*1000,2)}ms.")
print(f"Difference to z0: {jnp.linalg.norm(z0-moduli)}")
print(f"Difference to tau0: {jnp.abs(tau0-tau)}")

Beyond the hand-rolled loop, JAXVacua ships a JIT-compiled Newton solver as a method on FluxVacuaFinder — the subclass of FluxEFT that adds vacuum-finding utilities. The inner iteration runs inside jax.lax.while_loop under a @partial(jax.jit, static_argnums=...) boundary, so the full Newton refinement traces once and compiles to a single XLA kernel. The mode keyword selects which equation to solve — "SUSY" targets the \(F\)-term conditions \(D_I W = 0\) (relevant here), while the default mode=None solves the general extremum condition \(\partial_I V = 0\) used in the non-SUSY section below.

# FluxVacuaFinder inherits from FluxEFT; same fixture, same model_ID.
finder = jvc.FluxVacuaFinder(h12=h12, model_ID=1)

# Run the built-in JIT-compiled Newton solver in SUSY mode.
moduli_start = 2 * jnp.ones(h12) * 1j  # same starting point as the hand-rolled loop above
tau_start = 1j * 2.

moduli_n, tau_n, res = finder.newton_method_flux_vacua(
    moduli_start, tau_start, fluxes,
    mode="SUSY",
    tol=1e-12,
    max_iters=50,
    step_size_Newton=1.,
)
print(f"|D_I W| residual after JIT'd Newton: {res:.2e}")

The underlying _newton_method_flux_vacua_complex is decorated with @partial(jax.jit, static_argnums=(4, 5, 6, 7, 8)) — the static slots cover mode, step_size_Newton, tol, max_iters, and print_progress — so each unique combination of these triggers one trace+compile and is then cached. The cell above already paid the compile cost; the %%timeit below therefore reports the cached kernel’s runtime.

%%timeit
finder.newton_method_flux_vacua(
    moduli_start, tau_start, fluxes,
    mode="SUSY", tol=1e-12, max_iters=50,step_size_Newton=1.,
)

After jit compilation, the minimsation is very fast (of order a few milliseconds).

We can also employ alternative minimsation algorithms e.g. from the optimisation library scipy.optimize.. Here, we can use scipy.optimize.root as follows

x0  = jnp.array([0,1,0,1,0,1.])*2

result = root(model.DW_x,x0,args=(fluxes,), jac=model.dDW_x, method='hybr',tol=1e-10)

result

Let us check how close the solution is to the expected one, we compute

X = result.x
X = X[::2]+1j*X[1::2]

print(f"Difference to z0: {jnp.linalg.norm(z0-X[:2])}")
print(f"Difference to tau0: {jnp.abs(tau0-X[2])}")

We can now verify explicitly that the above solution corresponds to a minimum of the potential by looking at the Hessian:

z,cz,t,ct = model._convert_real_to_complex(result.x)
H = model.hessian(z,cz,t,ct,fluxes)
H_evals = jnp.linalg.eigvalsh(H)
H_evals

For a SUSY minimum, it can sometimes be convenient to use the faster implementation using mode="SUSY":

H_SUSY = model.hessian(z,cz,t,ct,fluxes,mode="SUSY")
H_SUSY_evals = jnp.linalg.eigvalsh(H_SUSY)
H_evals/H_SUSY_evals

We see that the eigenvalues are almost exactly the same. CAVEAT: Having said that, numerically setting \(D_I W=0\) in the code might lead to numerical errors that may become dangerous in siutations where high precision is required.

We can compare the time difference

%%timeit
H = jax.block_until_ready(model.hessian(z,cz,t,ct,fluxes))

%%timeit
H_SUSY = jax.block_until_ready(model.hessian(z,cz,t,ct,fluxes,mode="SUSY"))

We can compue the masses of the fields as follows

M = model.mass_matrix(z,cz,t,ct,fluxes)
masses = jnp.linalg.eigvalsh(M)
masses

We can again speed up the computation by using mode="SUSY"

M_SUSY = model.mass_matrix(z,cz,t,ct,fluxes,mode="SUSY")
masses_SUSY = jnp.linalg.eigvalsh(M_SUSY)
masses_SUSY/masses

This gives again the same output. We again find a significant speed-up:

%%timeit
M = jax.block_until_ready(model.mass_matrix(z,cz,t,ct,fluxes))

%%timeit
M_SUSY = jax.block_until_ready(model.mass_matrix(z,cz,t,ct,fluxes,mode="SUSY"))

Non-supersymmetric minima with \(D_IW\neq 0\)#

First, let us look at the following solution

# Moduli values
z0 = jnp.array([1.19 + 1.01j , 2.94 + 3.86j])

# Axio-dilaton values
tau0 = -0.40 + 5.30j

# Choices of fluxes
fluxes = jnp.array([0., 18., 3., 2., 4., -4., 0., -7., 0., 0., 0., 1.])

This solution was presented in appendix A of 2308.15525 (solution c in Tables 2 and 3). It corresponds to a dS minimum of the potential

model.dV(z0,jnp.conj(z0),tau0,jnp.conj(tau0),fluxes)

at which the \(F\)-terms are non-vanishing \(D_IW\neq 0\) as can be verified by computing

model.DW(z0,jnp.conj(z0),tau0,jnp.conj(tau0),fluxes)

Of course, the numerical tolerance \(|\partial_I V|<\epsilon\) is only as good as the precision used above which we improve with the numerical optimiser below.

We now select some starting point and try to find the above solution through numerical minimisation.

# Set starting point
moduli = 1+jnp.ones(h12)*1j
tau = 1j*1.

First, let us try Newton’s method

# Set parameters for minimsation
step_size_Newton = 0.1
tol = 1e-12

print("Jit compiling...")

# Get real representation 
x = model._convert_complex_to_real(moduli,jnp.conj(moduli),tau,jnp.conj(tau))
# Compute objective
dV = model.dV_x(x,fluxes)
# Compute Jacobian
ddV = model.ddV_x(x,fluxes)

print("Start minimisation...")
counter = 0
tic = time.time()
res = 10.
while res>tol:
    
    # Compute objective
    dV = model.dV_x(x,fluxes)
    # Compute Jacobian
    ddV = model.ddV_x(x,fluxes)
    # Compute step
    delta = -jnp.linalg.inv(ddV)@dV
    # Taking step in moduli and tau direction
    x = x+step_size_Newton*delta
    # Compute new residual
    res = jnp.sum(jnp.abs(dV))
    counter+=1
    
    print(f"residual: {res}                     counter: {counter}           ", end="\r")
    
    if counter>10**3:
        print("")
        print(f"STOPPED RUN: residual: {res}                     counter: {counter}           ")
        break
        
moduli,_,tau,_ = model._convert_real_to_complex(x)
toc = time.time()
print("")
print(f"Finished in {np.around((toc-tic)*1000,2)}ms.")
print(f"Difference to z0: {jnp.linalg.norm(z0-moduli)}")
print(f"Difference to tau0: {jnp.abs(tau0-tau)}")

The non-SUSY case calls the same solver — only mode changes. With mode=None (the default) the body inverts the Hessian of the scalar potential and follows dV instead of D_I W, so the iteration converges to a general extremum \(\partial_I V = 0\). Because mode is a static argument, this triggers a fresh JIT compile (a distinct cache entry from the SUSY call in §2.1).

# Same finder, same fluxes — only the starting point and `mode` differ.
moduli_start = 1+jnp.ones(h12) * 1j   # non-SUSY starting point used above
tau_start = 1j * 3.

moduli_n, tau_n, res = finder.newton_method_flux_vacua(
    moduli_start, tau_start, fluxes,
    mode=None,        # solve ∂_I V = 0
    tol=1e-12,
    max_iters=400,
    step_size_Newton=0.1,
    solver_mode="real"
)
print(f"|∂_I V| residual after JIT'd Newton: {res:.2e}")

%%timeit
finder.newton_method_flux_vacua(
    moduli_start, tau_start, fluxes,
    mode=None, tol=1e-12, max_iters=400,
    step_size_Newton=0.1, solver_mode="real"
)

After jit compilation, the minimsation is very fast (of order a few milliseconds).

We can also employ more sophisticated minimsation algorithms e.g. from the optimisation library scipy.optimize.. Here, we can use scipy.optimize.root as follows

eps=1e-10
x0=np.array([0,1,0,1,0,1.])*2

result = root(model.dV_x,x0,args=(fluxes,), jac=model.ddV_x, method='hybr',tol=eps)

result

Let us check how close the solution is to the expected one, we compute

X = result.x
X = X[::2]+1j*X[1::2]

print(f"Difference to z0: {jnp.linalg.norm(z0-X[:2])}")
print(f"Difference to tau0: {jnp.abs(tau0-X[2])}")

We can now verify explicitly that the above solution corresponds to a minimum of the potential by looking at the Hessian:

z,cz,t,ct = model._convert_real_to_complex(result.x)
H = model.hessian(z,cz,t,ct,fluxes)
H_evals = jnp.linalg.eigvalsh(H)
H_evals

We can compue the masses of the fields as follows

M = model.mass_matrix(z,cz,t,ct,fluxes)
masses = jnp.linalg.eigvalsh(M)
masses

Shared finder utilities#

The merged FluxVacuaFinder exposes a small set of stateless helpers that work on any converged candidate, regardless of which solver produced it. The most useful for post-hoc analysis is finder.classify_solution(x, flux), which returns a dict with the scalar potential value \(V\), \(|DW|\), the sorted Hessian eigenvalues, and SUSY/minimum flags. Below we apply it to the non-SUSY solution just computed.

# Rebuild the real-coord solution vector from the (moduli_n, tau_n) returned
# by Newton in the previous code cell.
x_sol = model._convert_complex_to_real_nondif(moduli_n,tau_n)

info = finder.classify_solution(x_sol, fluxes)
print(f"V          = {info['V']:+.3e}")
print(f"|DW|       = {info['|DW|']:.3e}")
print(f"is_susy    = {info['is_susy']}")
print(f"is_minimum = {info['is_minimum']}")
print(f"Nflux      = {info['Nflux']:.0f}")
print(f"eigvals    = {info['eigenvalues']}")

The returned dict has six keys:

V — scalar potential at the candidate.
|DW| — sum of \(|F\text{-terms}|\); flagged as SUSY when below \(10^{-6}\).
eigenvalues — sorted real Hessian eigenvalues.
is_susy — bool, |DW| < 1e-6.
is_minimum — bool, all eigenvalues exceed min_tol (default \(10^{-6}\)).
Nflux — \(|\mathrm{tadpole}(\mathrm{flux})|\), the D3-charge.

The same helper is also available as a module-level function jaxvacua.flux_utils.classify_solution(model, x, flux) for stateless callers. A batched, JIT-vmap’d variant finder._classify_batch(x_batch, flux_batch) is used internally by finder.sample_critical_points (see Notebook 07: ISD sampling) for fast classification of thousands of candidates in a single XLA call.

Projection on ISD and AISD fluxes#

The three-form flux \(G_3\) is naturally decomposed according to its Hodge type. It is well known that supersymmetric fluxes are of type \((2,1)\) and primitive, see 1707.01095 for references. In the low-energy effective theory, this condition is reflected by the fact that, at the vacuum point in complex structure moduli space, the flux vector \(\vec{N}=\vec{f}-\tau\vec{h}\) has non-vanishing components only along the \(D_i \Pi\) directions. Supersymmetry is already broken at tree level by the Kähler moduli directions by the presence of a non-vanishing \((0,3)\) piece (corresponding to \(W_0 \neq 0\)). The \((2,1)\)-components together with the \((0,3)\)-component furnish the ISD part of the complex 3-form \(G_3\).

In contrast, the non-supersymmetric solutions with \(D_IW\neq 0\) include, in addition to the \((2,1)\)-components, components along the \((1,2)\) and \((3,0)\) directions. They form the AISD part of the complex 3-form \(G_3\).

To project the integer fluxes in our symplectic basis onto the ISD and AISD components, recall that the set \(\{ \Omega, \bar{\Omega}, D_i \Omega, \bar{D}_{\bar{\imath}} \bar{\Omega} \}\) forms a basis of three-forms for generic values of the moduli, see Candelas, de la Ossa 1990. The corresponding set \(\{ \Pi, \Pi^*, D_i \Pi, \bar{D}_{\bar{\imath}} \Pi^* \}\) constitutes an orthogonal basis with respect to the symplectic inner product.As a result, any vector \(\vec{V}\) in the period/charge vector space can be expanded as (see 1707.01095):

\[ \vec{V} = \mathrm{i} e^{K_{\text{cs}}} \Bigl [(\Pi^\dagger \Sigma \vec{V}) \Pi- (\Pi^T \Sigma \vec{V}) \Pi^*- (K^{\bar{\jmath} l} \bar{D}_{\bar{\jmath}} \Pi^\dagger \Sigma \vec{V}) D_l \Pi - (K^{j \bar{l}} D_j \Pi^T \Sigma \vec{V}) \bar{D}_{\bar{l}} \Pi^* \Bigl ] \]

where \(\Pi\) is evaluated at a generic point in complex structure moduli space. For the flux vector \(\vec{N}\), the expansion coefficients are given by

\[ \vec{N} = \mathrm{i} e^{K_{\text{cs}}} \left(\frac{D_{\bar{\tau}}\overline{W}}{K_{\bar{\tau}}} \Pi- W\Pi^*- \frac{(D_{\bar{\jmath}}D_{\bar{\tau}}\overline{W})}{K_{\bar{\tau}}}\, K^{i\bar{\jmath}} \, D_i \Pi- (D_i W) K^{i\bar{\jmath}} \bar{D}_{\bar{\jmath}} \Pi^*\right) \]

Let us compute the corresponding projections for the supersymmetric minimum from above:

# Moduli values
z_SUSY = jnp.array([0.3168651582886485 +3.2301486028589927j,
                    0.37113099110064307+3.0746240438915935j])

# Axio-dilaton values
tau_SUSY = -0.29834960566944685+6.803576151333438j

# Choices of fluxes
fluxes_SUSY = jnp.array([4.,  -3.,  -2.,  -2.,   3.,   2.,  39., -13.,  -4.,   0.,  -0.,  0.])

N_3_0,N_2_1,N_1_2,N_0_3 = model.projection_fluxes(z_SUSY,tau_SUSY,fluxes_SUSY,mode=None)

We can check that the \((3,0)\)- and \((1,2)\)-components vanish

np.sum(np.abs(N_3_0)),np.sum(np.abs(N_1_2))

The remaining components should therefor add up to the flux vectors fluxes_SUSY from above

np.around(N_2_1+N_0_3,10)

or more explicitly

np.sum(np.abs(N_2_1+N_0_3-fluxes_SUSY))

Next, let us consider the non-SUSY minimum from above

# Moduli values
z_nonSUSY = jnp.array([1.190671013458699 +1.0138094643117892j,2.9377413761449915+3.8593967759445964j])

# Axio-dilaton values
tau_nonSUSY = -0.403614509914007+5.301822679470445j

# Choices of fluxes
fluxes_nonSUSY = jnp.array([0., 18., 3., 2., 4., -4., 0., -7., 0., 0., 0., 1.])

We compute the projections as above

N_3_0,N_2_1,N_1_2,N_0_3 = model.projection_fluxes(z_nonSUSY,tau_nonSUSY,fluxes_nonSUSY,mode=None)

This time, the \((3,0)\)- and \((1,2)\)-components are non-vanishing

np.sum(np.abs(N_3_0)),np.sum(np.abs(N_1_2)),np.sum(np.abs(N_2_1)),np.sum(np.abs(N_0_3))

We note however that the \((3,0)\)- and \((1,2)\)-components are comparatively small to the ISD components. This is because the SUSY breaking scale is rather small (see 2308.15525 for a discussion) as can be seen from the value of the scalar potential at the minimum:

model.V(z_nonSUSY,jnp.conj(z_nonSUSY),tau_nonSUSY,jnp.conj(tau_nonSUSY),fluxes_nonSUSY)

We can again make the consistency check that the components add up to the integer quantised flux from above

np.around(N_3_0+N_2_1+N_1_2+N_0_3,10)