Theories and solvers

Note

Throughout, hatted symbols indicate quantities expressed in in-code dimensionless units: e.g. \(\hat{\phi}\) indicates the rescaled UV field \(\phi/M_{f1}\), where \(M_{f1}\) is some mass.

Base class

class phienics.solver.Solver(fem, source, fields, Mn=None, abs_du_tol=1e-08, rel_du_tol=1e-08, abs_res_tol=1e-10, rel_res_tol=1e-10, max_iter=100, criterion='residual', norm_change='linf', norm_res='linf')

Base class for solvers of models of screening.

This solver implements the Newton method given the equation for the iterations: if \(\mathcal{F}(u)=0\) is the original nonlinear equation (or system of equations), this solver iterates:

(1)\[0 = \mathcal{F}[u^{(k)}] + \int d{\bf x} \frac{\delta \mathcal{F}}{\delta u({\bf x})}[u^{(k)}] \left(u_{\rm true}-u^{(k)}\right)\]

at every iteration \(k\), until convergence. Eq. (1) must be supplied as linear_solver(). For further details see the paper .

Two convergence tests are available:

1. ‘residual’ (default):

(2)\[|| \mathcal{F}[u^{(k)}] || \leq \epsilon_{\rm rel}^{\rm (R)} || \mathcal{F}[u^{(0)}] || + \epsilon_{\rm abs}^{\rm (R)}\]

2. ‘change’, i.e. change in the solution:

(3)\[|| u^{(k)}-u^{(k-1)} || \leq \epsilon_{\rm rel}^{\rm (S)} || u^{(0)} || + \epsilon_{\rm abs}^{\rm (S)}\]

where \(\epsilon_{\rm rel}^{\rm (\cdot)}\) and \(\epsilon_{\rm abs}^{\rm (\cdot)}\) are relative and absolute tolerances, \(u^{(0)}\) is the initial guess, and \(|| \cdot ||\) is some norm on the solution space.

Arguments

fem

a fem.Fem instance

source

a source.Source instance

fields

a UV.UVFields or IR.IRFields instance

Mn

distances are rescaled following \(\hat{r} = M_n r\) inside the code. Default: \(M_n={r_S}^{-1}\)

abs_du_tol

if criterion is ‘change’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (S)}\) in Eq. (3), i.e. the absolute tolerance on the change in the solution between two iterations

rel_du_tol

if criterion is ‘change’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (S)}\) in Eq. (3), i.e. the relative tolerance on the change in the solution between two iterations

abs_res_tol

if criterion is ‘residual’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (R)}\) in Eq. (2), i.e. the absolute tolerance on the residuals at a given iteration

rel_res_tol

if criterion is ‘residual’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (R)}\) in Eq. (2), i.e. the relative tolerance on the residuals at a given iteration

max_iter

number of maximum iterations allowed. If the solver does not converge in max_iter iterations, it will stop and print a warning; however, it will still return the output and it will not throw an exception

criterion

‘change’ or ‘residual’ (default: ‘residual’): convergence criterion

norm_change

norm used in the evaluation of the ‘change’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

norm_res

norm used in the evaluation of the ‘residual’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

On(n, method='derivative', rescale=1.0, output_rescaled_op=False)

Computes the operators:

\[O_n = (-1)^{n+1} \frac{2n+2}{2n+1} \binom{3n}{n} \alpha^{2n+1} \lambda^n \nabla^2 \left((\nabla^2 \varphi)^{2n+1}\right) M^{-6n-2}\]

associated to a field \(\varphi\).

In the examples provided within \(\varphi\mathrm{enics}\), \(\varphi\) is the the UV field \(\phi\) or the IR field \(\pi\).

The operators are computed starting from the field’s rescaled Laplacian \(y\equiv\hat{\nabla}^2\hat{\varphi}\), so the field’s equation of motion must have been solved prior to invoking this method.

For the method used to compute \(\nabla^2 \left((\nabla^2 \varphi)^{2n+1}\right)\), see the documentation of Qn().

Arguments

n

order of the operator \(O_n\)

method

‘derivative’, ‘auxiliary’ or ‘gradient’ - see func:Qn

rescale

(optional) temporary auxiliary variable used to rescale the Laplacian during the computation of \(\nabla^2 \left((\nabla^2 \varphi)^{2n+1}\right)\), to prevent hitting the maximum/minimum representable number

output_rescaled_op

True to obtain the rescaled operator \(\hat{O}_n \equiv \hat{\nabla}^2 \left((\hat{\nabla}^2 \hat{\varphi})^{2n+1}\right)\) alongside the physical operator \(O_n\), in a tuple \((\hat{O}_n, O_n)\); False otherwise. Default: False.

Qn(n, method='derivative', rescale=1.0, output_rescaled_op=False)

Computes the operators:

\[Q_n = \frac{\alpha^n}{M^{3n-1}} \nabla^2( ( \nabla^2\varphi )^n )\]

associated to a field \(\varphi\).

In the examples provided within \(\varphi\mathrm{enics}\), \(\varphi\) is the the UV field \(\phi\) or the IR field \(\pi\).

The key to computing \(Q_n\) is obtaining the rescaled operator \(\hat{Q}_n \equiv \hat{\nabla}^2 \left((\hat{\nabla}^2 \hat{\varphi})^{n}\right)\). The Laplacian \(y\equiv\hat{\nabla}^2\hat{\varphi}\) is obtained from the solution to the equation of motion; \(\hat{Q}_n\) can then be computed in three different ways, specified by the method option:

• derivative (default):

  the rescaled operator is decomposed as:

  \[\hat{\nabla}^2( y^{n} ) = n y^{n-1} \hat{\nabla}^2 y + n (n-1) y^{n-2} \hat{\nabla} y \cdot \hat{\nabla} y\]

  before being projected on the discontinuous function space specified in the fem instance;

• auxiliary:

  \(w \equiv y^n\) is projected onto the discontinuous function space specified in the fem instance, then the code solves for the linear system \(\hat{Q}_n = \nabla^2(w)\);

• gradient:

  for the UV and IR theories supplied as \(\varphi\mathrm{enics}\) examples, the weak formulation of the equation \(\hat{Q}_n = \hat{\nabla}^2( y^n )\) is

  \[\int \hat{Q}_n v \hat{r}^2 d\hat{r} = \int \hat{\nabla}^2 ( y^n ) v \hat{r}^2 d\hat{r} = - \int \hat{\nabla}( y^n ) \hat{\nabla} v \hat{r}^2 d\hat{r} = - n \int y^{n-1} \hat{\nabla}y \hat{\nabla}v \hat{r}^2 d\hat{r}\]

  where \(v\) is a test function, so that the \(\hat{Q}_n\) operator can be obtained by solving

  \[\int \hat{Q}_n v \hat{r}^2 d\hat{r} = - n \int y^{n-1} \hat{\nabla}y \hat{\nabla}v \hat{r}^2 d\hat{r}\]

  which is the case for the gradient option.

Further details are given in the paper .

Note

The \(Q_n\) operators are computed starting from the field’s rescaled Laplacian \(y\equiv\hat{\nabla}^2\hat{\varphi}\), so the method solve() (solving the field’s equation of motion) must have been called before calling this method.

Arguments

n

order of the operator \(Q_n\)

method

‘derivative’, ‘auxiliary’ or ‘gradient’

rescale

(optional) temporary auxiliary rescaling for \(Q_n\), used in tests or to prevent hitting the maximum/minimum representable number

output_rescaled_op

True to obtain the rescaled operator \(\hat{Q}_n \equiv \hat{\nabla}^2 \left((\hat{\nabla}^2 \hat{\varphi})^{2n+1}\right)\) alongside the physical operator \(Q_n\), in a tuple \((\hat{O}_n, O_n)\); False otherwise. Default: False.

compute_errors(du_k, u_k)

Computes measures of error at a given iteration.

With reference to Eq. (2) and Eq. (3), this function computes the norm of the residual \(|| \mathcal{F}[u^{(k)}] ||\) and the norm of the change in the solution \(|| u^{(k)}-u^{(k-1)} ||\) at one given iteration.

Parameters

du_k

difference in the solution over the whole domain at this iteration.

u_k

solution at this iteration, for the computation of the residual \(||F(u_k)||\)

grad(field, radial_units)

Returns the gradient of the scalar field in input.

The gradient is computed by projecting the radial derivative of the field onto the discontinuous function space specified in solver.fem.

Arguments

field

the field of which you want to compute the gradient

radial_units

‘physical’ or ‘rescaled’: the former returns the field’s gradient with respect to physical distances (units \({M_p}^{-1}\)), the latter returns the gradient with respect to the dimensionless rescaled distances used within the code

scalar_force(field)

Returns the magnitude of the scalar force associated to the input field, per unit mass (units \(M_p\)):

\[F_{\varphi} = \frac{\nabla\varphi}{M_P}\]

if \(\varphi\) is the input field.

Arguments field

the field associated to the scalar force

solve()

Solves the non-linear equation iteratively using the Newton’s method, starting from a linear solver as in Eq. (1).

strong_residual(sol, units='rescaled', norm='linf')

Returns the residual, either as a function on the domain or as a norm of that function.

The former is just a wrapper around strong_residual_form(), which is determined by the model (see e.g. UV.UVSolver.strong_residual_form()); the latter computes a choice of norm defined by the parameter norm.

Arguments

sol

the solution \(u\) on which the residual \(\mathcal{F}[u]\) is to be computed

units

‘rescaled’ (for the rescaled units used inside the code) or ‘physical’, for physical units

norm

‘L2’, ‘linf’ or ‘none’. If ‘L2’ or ‘linf’: compute the \(L_2\) or \(\ell_{\infty}\) norm of the residual; if ‘none’, return the full function - as opposed to its norm.

UV theory

class phienics.UV.UVFields(m=1e-51, M=1e-48, Mp=1.0, alpha=0.4, lam=0.7)

Physical parameters of the UV-complete theory described in Section VI of De Rham et al:

\[ \begin{align}\begin{aligned}&\Box\phi - m^2\phi - \alpha\Box H = \frac{\rho}{M_p}\\&\Box H - M^2 H - \alpha\Box\phi - \frac{\lambda}{3!} H^3 = 0\end{aligned}\end{align} \]

plus an additional coupling to matter to study screening.

\(\phi\) is a cosmological light scalar field of mass \(m\), \(H\) is a self-interacting more massive scalar field of mass \(M\), \(\rho\) is a massive source to which \(\phi\) is coupled, \(\alpha\) is the coupling between the two scalar fields, \(\lambda\) is the strength of the self-coupling of \(H\) and \(M_p\) is the Planck mass.

Arguments

m

mass of the light field \(m\)

M

mass of the heavy field \(M\)

Mp

Planck mass \(M_p\)

alpha

coupling constant \(\alpha\)

lam

self-coupling \(\lambda\) (note that ‘lambda’ is a reserved word in Python, and can’t be used as a variable name)

class phienics.UV.UVSolver(fem, source, fields, Mn=None, Mf1=None, Mf2=None, abs_du_tol=1e-08, rel_du_tol=1e-08, abs_res_tol=1e-10, rel_res_tol=1e-10, max_iter=100, criterion='residual', norm_change='linf', norm_res='linf')

Solver for the UV-complete theory described in Section VI of De Rham et al (plus a coupling to matter). Derives from solver.Solver.

For better accuracy in the computation of the fields’ Laplacians, within the code the original system of equations is expressed as:

\[\begin{split}& Y - m^2 \phi - \alpha Z = \frac{\rho}{M_P} \\ & Z - M^2 H - \alpha Y - \frac{\lambda}{6} H^3 = 0 \\ & \nabla^2 \phi = Y \\ & \nabla^2 H = Z \end{split}\]

Arguments

fem

a fem.Fem instance

source

a source.Source instance

fields

a UVFields instance

Mn

distances are rescaled following \(\hat{r} = M_n r\) inside the code. Default: \(M_n={r_S}^{-1}\)

Mf1

field rescaling for \(\phi\): \(\hat{\phi}=\phi/M_{f1}\). Default: \(M_s M_n / M_P\).

Mf2

field rescaling for \(H\): \(\hat{H}=H/M_{f2}\). Default: \(|\alpha| M_{f1}\).

abs_du_tol

if criterion is ‘change’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (S)}\) in Eq. (3), i.e. the absolute tolerance on the change in the solution between two iterations

rel_du_tol

if criterion is ‘change’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (S)}\) in Eq. (3), i.e. the relative tolerance on the change in the solution between two iterations

abs_res_tol

if criterion is ‘residual’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (R)}\) in Eq. (2), i.e. the absolute tolerance on the residuals at a given iteration

rel_res_tol

if criterion is ‘residual’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (R)}\) in Eq. (2), i.e. the relative tolerance on the residuals at a given iteration

max_iter

number of maximum iterations allowed. If the solver does not converge in max_iter iterations, it will stop and print a warning; however, it will still return the output and it will not throw an exception

criterion

‘change’ or ‘residual’ (default: ‘residual’): convergence criterion

norm_change

norm used in the evaluation of the ‘change’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

norm_res

norm used in the evaluation of the ‘residual’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

compute_physics(sol)

Computes field profiles and gradients in physical units, starting from the dimensionless profiles that are obtained solving the equations of motion with solve. Then computes the scalar force associated to \(\phi\) (physical units only).

get_Dirichlet_bc()

Returns the Dirichlet boundary conditions for the UV theory:

\[\hat{\phi}(\hat{r}_{\rm max}) = \hat{H}(\hat{r}_{\rm max}) = \hat{Y}(\hat{r}_{\rm max}) = \hat{Z}(\hat{r}_{\rm max}) = 0\]

The Neumann boundary conditions

\[\frac{d\hat{\phi}}{d\hat{r}}(0) = \frac{d\hat{H}}{d\hat{r}}(0) = 0\]

are natural boundary conditions in the finite element method, and are implemented within the variational formulation.

initial_guess()

Obtains an initial guess for the Newton solver.

This is done by solving the equation of motion for \(\lambda=0\), which form a linear system. All other parameters are unchanged.

linear_solver(u_k)

Solves the (linear) Newton iteration Eq. (1) for the UV theory.

For the UV theory, the form of the Newton iteration is:

\[ \begin{align}\begin{aligned}& \int \hat{Y} v_1 \hat{r}^2 d\hat{r} - \left(\frac{m}{M_n}\right)^2 \int \hat{\phi} v_1 \hat{r}^2 d\hat{r} - \alpha \int \hat{Z} v_1 \hat{r}^2 d\hat{r} + \\& + \int \hat{Z} v_2 \hat{r}^2 d\hat{r} - \left(\frac{M}{M_n}\right)^2 \int \hat{H} v_2 \hat{r}^2 d\hat{r} - \alpha \int \hat{Y} v_2 \hat{r}^2 d\hat{r} - \frac{\lambda}{2} \int {\hat{H}_k}^2 \hat{H} v_2 \hat{r}^2 d\hat{r} +\\& - \int \hat{\nabla}\hat{\phi} \cdot \hat{\nabla} v_3 \hat{r}^2 d\hat{r} - \int \hat{Y} v_3 \hat{r}^2 d\hat{r} - \int \hat{\nabla}\hat{H} \cdot \hat{\nabla} v_4 \hat{r}^2 d\hat{r} - \int \hat{Z} v_4 \hat{r}^2 d\hat{r} = \\& = \frac{M_n}{M_{f1}} \int \frac{\hat{\rho}}{M_P} v_1 \hat{r}^2 d\hat{r} - \int \frac{\lambda}{3} {\hat{H}_k}^3 v_2 \hat{r}^2 d\hat{r}\end{aligned}\end{align} \]

Arguments

u_k

solution at the previous iteration

output_term(eqn=1, term='LHS', norm='none', units='rescaled', output_label=False)

Outputs the different terms in the UV system of equations.

The terms can be output in physical or dimensionless in-code units, by choosing either units=’physical’ or units=’rescaled’.

Possible choices of terms are listed in the tables below; recall that \(Y\equiv\nabla^2\phi\) and \(Z\equiv\nabla^2 H\).

• eqn=1:

units=	term=’LHS’	term=’RHS’
‘rescaled’	\(\hat{Y} - \left(\frac{m}{M_n}\right)^2\hat{\phi} - \frac{M_{f2}}{M_{f1}}\alpha\hat{Z}\)	\(\frac{\hat{\rho}}{M_P}\frac{Mn}{M_{f1}}\)
‘physical’	\(Y -m^2 \phi -\alpha Z\)	\(\frac{\rho}{M_P}\)

units=	term=1	term=2	term=3	term=4
‘rescaled’	\(\hat{Y}\)	\(-\left(\frac{m}{M_n}\right)^2\hat{\phi}\)	\(-\frac{M_{f2}}{M_{f1}}\alpha\hat{Z}\)	same as ‘RHS’
‘physical’	\(Y\)	\(-m^2\phi\)	\(-\alpha Z\)	same as ‘RHS’

• eqn=2:

units=	term=’LHS’	term=’RHS’
‘rescaled’	\(\hat{Z} -\left(\frac{M}{M_n}\right)^2\hat{H} - \alpha\frac{M_{f1}}{M_{f2}}\hat{Y}\)	\(-\frac{\lambda}{6} \left( \frac{M_{f2}}{M_n} \right)^2 \hat{H}^3\)
‘physical’	\(Z -M^2 H -\alpha Y\)	\(-\frac{\lambda}{6} H^3\)

units=	term=1	term=2	term=3	term=4
‘rescaled’	\(\hat{Z}\)	\(-\left(\frac{M}{Mn}\right)^2\hat{H}\)	\(-\alpha\frac{M_{f1}}{M_{f2}}\hat{Y}\)	same as ‘RHS’
‘physical’	\(Z\)	\(-M^2 H\)	\(-\alpha Y\)	same as ‘RHS’

Although the right-hand-side (RHS) of the equation of motion for \(H\) is 0, this choice allows better comparison of the scale of the equation against the residuals.

• eqn=3:

units=	term=’LHS’	term=1	term=2
‘rescaled’	\(\hat{\nabla}^2\hat{\phi}-\hat{Y}\)	\(\hat{\nabla}^2\hat{\phi}\)	\(\hat{Y}\)
‘physical’	\(\nabla^2\phi - Y\)	\(\nabla^2\phi\)	\(Y\)

• eqn=4:

units=	term=’LHS’	term=1	term=2
‘rescaled’	\(\hat{\nabla}^2\hat{H}-\hat{Z}\)	\(\hat{\nabla}^2\hat{H}\)	\(\hat{Z}\)
‘physical’	\(\nabla^2 H - Z\)	\(\nabla^2 H\)	\(Z\)

Equations (1) and (2) allow to look at the terms in the equations of motion, Eq. (3) and (4) enforce consistency relations \(Y = '\nabla^2\phi\) and \(Z = \nabla^2 H\).

Note

In this method, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Parameters

eqn

(integer) choice of equation

term

choice of LHS, RHS or a number; for eqn=1,2, valid terms range from 1 to 4, for eqn=3,4, valid terms range from 1 to 2

norm

‘L2’, ‘linf’ or ‘none’. If ‘L2’ or ‘linf’: compute the \(L_2\) or \(\ell_{\infty}\) norm of the residual; if ‘none’, return the full term over the box - as opposed to its norm.

units

rescaled (default) or physical; choice of units for the output

output_label

if True, output a string with a label for the term (which can be used, e.g. in plot legends)

strong_residual_form(sol, units)

Computes the residual with respect to the strong form of the equations.

The total residual is obtained by summing the residuals of all equations:

\[F = F_1 + F_2 + F_3 + F_4\]

where, in dimensionless in-code units (units=’rescaled’):

\[ \begin{align}\begin{aligned}& F_1(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \hat{Y} - \left( \frac{m}{M_n} \right)^2\hat{\phi} - \alpha \frac{M_{f2}}{M_{f1}} \hat{Z} - \frac{\hat{\rho}}{M_p}\frac{M_n}{M_{f1}}\\& F_2(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \hat{Z} - \left( \frac{M}{M_n} \right)^2 \hat{H} - \alpha \frac{M_{f1}}{M_{f2}} \hat{Y} - \frac{\lambda}{6} \left( \frac{M_{f2}}{M_n} \right)^2 \hat{H}^3\\& F_3(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \hat{\nabla}^2\hat{\phi} - \hat{Y}\\& F_4(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \hat{\nabla}^2\hat{H} - \hat{Z}\end{aligned}\end{align} \]

and in physical units (units=’physical’):

\[ \begin{align}\begin{aligned}& F_1(\phi,H,Y,Z) = Y - m^2 \phi -\alpha H - \frac{\rho}{M_P}\\& F_2(\phi,H,Y,Z) = Z - M^2 H - \alpha\phi - \frac{\lambda}{6} H^3\\& F_3(\phi,H,Y,Z) = \nabla^2\phi - Y\\& F_4(\phi,H,Y,Z) = \nabla^2 H - Z\end{aligned}\end{align} \]

Note

In this function, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Note that the weak residual in weak_residual_form() is just the scalar product of the strong residuals by test functions.

Parameters

sol

the solution with respect to which the weak residual is computed.

units

‘rescaled’ (for the rescaled units used inside the code) or ‘physical’, for physical units

weak_residual_form(sol)

Computes the residual with respect to the weak form of the equations:

\[F = F_1 + F_2 + F_3 + F_4\]

with

\[ \begin{align}\begin{aligned}\begin{split}&F_1(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \int\hat{Y} v_1 \hat{r}^2 d\hat{r} - \left(\frac{m}{M_n}\right)^2 \int\hat{\phi} v_1 \hat{r}^2 d\hat{r} - \alpha \frac{M_{f2}}{M_{f1}} \int \hat{Z} v_1 \hat{r}^2 d\hat{r} - \frac{M_n}{M_{f1}} \int \frac{\hat{\rho}}{M_p} v_1 \hat{r}^2 d\hat{r} \\\end{split}\\\begin{split}&F_2(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = \int \hat{Z} v_2 \hat{r}^2 d\hat{r} - \left( \frac{M}{M_n} \right)^2 \int \hat{H} v_2 \hat{r}^2 d\hat{r} - \alpha \frac{M_{f1}}{M_{f2}} \int \hat{Y} v_2 \hat{r}^2 d\hat{r} - \frac{\lambda}{6} \left( \frac{M_{f2}}{M_n} \right)^2 \int \hat{H}^3 v_2 \hat{r}^2 d\hat{r} \\\end{split}\\\begin{split}&F_3(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = - \int \hat{\nabla} \hat{\phi} \cdot \hat{\nabla} v_3 \hat{r}^2 d\hat{r} - \hat{Y} v_3 \hat{r}^2 d\hat{r} \\\end{split}\\&F_4(\hat{\phi},\hat{H},\hat{Y},\hat{Z}) = - \int \hat{\nabla} \hat{H} \cdot \hat{\nabla} v_4 \hat{r}^2 d\hat{r} - \hat{Z} v_4 \hat{r}^2 d\hat{r}\end{aligned}\end{align} \]

The weak residual is employed within solver.Solver.solve() to check convergence - also see solver.Solver.compute_errors().

Parameters

sol

the solution with respect to which the weak residual is computed.

IR theory

class phienics.IR.IRFields(m=1e-51, Lambda=1e-49, Mp=1.0, epsilon=0.02, n=3)

Physical parameters of the massive Galileon theory of equation of motion:

(4)\[\Box\pi -m^2\pi - \frac{\epsilon}{\Lambda^{3n-1}} \Box\left( \Box\pi^n \right) = \frac{\rho}{M_P}\]

Here \(\pi\) is a massive scalar field of mass \(m\), \(\Lambda\) is a cut-off scale with units of mass, and \(\rho\) is the energy density of the source. The dimensionless parameter \(\epsilon\) regulates the strength of the nonlinear term which is expected to give rise to Vainshtein screening.

Arguments

m

mass of Galileon field \(\pi\)

Lambda

cut-off scale of the theory

Mp

Planck mass \(M_P\)

epsilon

nonlinear parameter \(\epsilon\)

n

exponent \(n\) in the theory

class phienics.IR.IRSolver(fem, source, fields, Mn=None, Mf1='NL', guess_choice='NL', abs_du_tol=1e-08, rel_du_tol=1e-08, abs_res_tol=1e-10, rel_res_tol=1e-10, max_iter=100, criterion='residual', norm_change='linf', norm_res='linf')

Solver for the Galileon theory in Eq. (4). Derives from solver.Solver.

For better accuracy in the computation of the field’s Laplacian and of the nonlinear term, the equation of motion Eq. (4) is expressed within the code as:

\[\begin{split}& \nabla^2 \phi = Y \\ & W = Y^n \\ & Y - m^2 \phi - \frac{\epsilon}{\Lambda^{3n-1}} \nabla^2 W = \frac{\rho}{M_P}\end{split}\]

Two off-the-shelf choices are available for the initial guess:

• guess_choice=’NL’ (default):

  this choice assumes that the nonlinear term in the equation of motion is dominant, obtaining the approximation

  \[-\nabla^2 W_0 \approx \frac{\Lambda^{3n-1}}{\epsilon} \frac{\rho}{M_P}\]

  from which \(Y\) is obtained as \(Y_0=\sqrt[n]{W_0}\), and finally \(\nabla^2\pi_0=Y_0\)

• guess_choice=’KG’:

  this choice assumes that the nonlinear term in the equation of motion is subdominant, obtaning the approximation

  \[\nabla^2\pi - m^2 \pi \approx \frac{\rho}{M_P}\]

  which is a linear system and can be solved as is.

Two off-the-shelf choice of dimensionless units are similarly available for the scalar field \(\pi\), when the user does not wish to set it manually; having defined \(\hat{\pi}=\pi/M_{f1}\), one can choose:

• Mf1=’NL’ (default):

  with this choice, the code attempts to find an optimal rescaling for \(\pi\) based on the assumption that the nonlinear term is dominant, by setting

  \[M_{f1} = \left( \Lambda / M_n \right)^{(3n-1)/n} ( M_S / M_P)^{1/n} \, m\]

• Mf1=’source’:

  this choice results in a dimensionless source term that is \(\approx O(1)\), by setting:

  \[M_{f1} = M_S / M_P M_n\]

where \(M_S\) is the source mass and, in both cases, \(M_n\) is a unit of inverse distance used in the definition of the in-code dimensionless radial distances: \(\hat{r}\equiv r M_n\) (for the definition of other symbols, please see IRFields).

Arguments

fem

a fem.Fem instance

source

a source.Source instance

fields

a UVFields instance

Mn

distances are rescaled following \(\hat{r} = M_n r\) inside the code. Default: \(M_n={r_S}^{-1}\)

Mf1

field rescaling for \(\pi\): \(\hat{\pi}=\pi/M_{f1}\); choice of ‘NL’, ‘source’ or a real number. Default:Mf1=’NL’

guess_choice

choice of NL or KG. Default: NL

abs_du_tol

if criterion is ‘change’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (S)}\) in Eq. (3), i.e. the absolute tolerance on the change in the solution between two iterations

rel_du_tol

if criterion is ‘change’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (S)}\) in Eq. (3), i.e. the relative tolerance on the change in the solution between two iterations

abs_res_tol

if criterion is ‘residual’, abs_du_tol is the quantity \(\epsilon_{\rm abs}^{\rm (R)}\) in Eq. (2), i.e. the absolute tolerance on the residuals at a given iteration

rel_res_tol

if criterion is ‘residual’, rel_du_tol is the quantity \(\epsilon_{\rm rel}^{\rm (R)}\) in Eq. (2), i.e. the relative tolerance on the residuals at a given iteration

max_iter

number of maximum iterations allowed. If the solver does not converge in max_iter iterations, it will stop and print a warning; however, it will still return the output and it will not throw an exception

criterion

‘change’ or ‘residual’ (default: ‘residual’): convergence criterion

norm_change

norm used in the evaluation of the ‘change’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

norm_res

norm used in the evaluation of the ‘residual’ criterion (default: ‘linf’, i.e. maximum absolute value at nodes)

KG_initial_guess()

Obtains an initial guess for the Galileon equation of motion by assuming the nonlinear term is subdominant, i.e.:

\[\Box\pi - m^2\pi \approx \frac{\rho}{Mp}\]

The initial guess is computed by first solving the system of equations:

\[\begin{split}& \hat{\nabla}^2\hat{\pi} = \hat{Y} \\ & \hat{Y} - \left( \frac{m}{M_n} \right)^2\pi = \frac{\hat{\rho}}{M_p}\end{split}\]

and then obtaining \(\hat{W}=\hat{Y}^n\) by projection.

NL_initial_guess(y_method='vector')

Obtains an initial guess for the Galileon equation of motion by assuming the nonlinear term is dominant, i.e.:

\[-\frac{\epsilon}{\Lambda^{3n-1}}\nabla^2(\nabla^2\pi^n) \approx \frac{\rho}{M_P}\]

The initial guess is computed by first solving the Poisson equation:

\[-\hat{\nabla}\hat{W} =\left( \frac{\Lambda}{M_n} \right)^{3n-1} \left(\frac{M_n}{M_{f1}}\right)^n \frac{\hat{\rho}}{\epsilon M_P}\]

and then obtaining \(\hat{Y}=\sqrt[n]{\hat{W}}\) through one of three methods explained below. Finally, \(\hat{\pi}\) is computed by solving the Poisson equation

\[\hat{\nabla}\hat{\pi} = \hat{Y}.\]

The main methods to obtain \(\hat{Y}\) from \(\hat{Z}\) are interpolation and projection: the standard FEniCS implementation for both can be chosen by setting y_method=’interpolate’ and y_method=’project’.

A third method (y_method=’vector’, default), formally identical to interpolation, consists in assigning \(\hat{Y}\)’s value at all nodes through e.g.:

y.vector().set_local( np.sqrt( np.abs( w.vector().get_local() ) ) )

However, because of differences in the implementations of \(\sqrt[n]{\cdot}\) called by the two methods, the latter generally gives better results compared to ‘interpolate’.

Arguments

y_method

‘vector’ (default), ‘interpolate’ or ‘project’

get_Dirichlet_bc()

Returns the Dirichlet boundary conditions for the IR theory:

\[\hat{\pi}(\hat{r}_{\rm max}) = \hat{W}(\hat{r}_{\rm max}) = \hat{Y}(\hat{r}_{\rm max}) = 0\]

The Neumann boundary conditions

\[\left\lbrace \frac{d\hat{\pi}}{d\hat{r}}(0) = 0; \frac{d\hat{W}}{d\hat{r}}(0) = \mathrm{finite} \right\rbrace\]

are natural boundary conditions in the finite element method, and are implemented within the variational formulation.

linear_solver(u_k)

Solves the (linear) Newton iteration Eq. (1) for the IR theory.

For the IR theory, the form of the Newton iterations is:

\[ \begin{align}\begin{aligned}& - \int \hat{\nabla}\hat{\pi} \cdot \hat{\nabla} v_1 \hat{r}^2 d\hat{r} - \int \hat{Y} v_1 \hat{r}^2 d\hat{r} +\\& + \int \hat{W} v_2 \hat{r}^2 d\hat{r} -n \int {\hat{Y}_k}^{n-1} \hat{Y} v_2 \hat{r}^2 d\hat{r} +\\& + \int \hat{Y} v_3 \hat{r}^2 d\hat{r} - \int \left( \frac{m}{M_n} \right)^2 \hat{\pi} v_3 \hat{r}^2 d\hat{r} + \epsilon \left( \frac{M_n}{\Lambda} \right)^{3n-1} \left( \frac{M_{f1}}{M_n} \right)^{n-1} \int \nabla\hat{W} \cdot \nabla v_3 \hat{r}^2 d\hat{r}\\& = (1-n) \int {\hat{Y}_k}^n v_2 \hat{r}^2 d\hat{r} + \int \frac{\hat{\rho}}{M_P} \frac{M_n}{M_{f1}} v_3 \hat{r}^2 d\hat{r}\end{aligned}\end{align} \]

Arguments

u_k

solution at the previous iteration

output_term(eqn=1, term='LHS', norm='none', units='rescaled', output_label=False)

Outputs the different terms in the IR system of equations.

The terms can be output in physical or dimensionless in-code units, by choosing either units=’physical’ or units=’rescaled’.

Possible choices of terms are listed in the tables below; recall that \(Y\equiv\nabla^2\phi\) and \(W \equiv Y^n\).

• eqn=1:

units=	term=’LHS’	term=1	term=2
‘rescaled’	\(\hat{Y} - \hat{\nabla}^2\hat{\pi}\)	\(\hat{\nabla}^2\hat{\pi}\)	\(\hat{Y}\)
‘physical’	\(Y - \nabla^2\pi\)	\(\nabla^2\pi\)	\(Y\)

• eqn=2:

units=	term=’LHS’	term=1	term=2
‘rescaled’	\(\hat{W} - \hat{Y}^n\)	\(\hat{Y}^n\)	\(\hat{W}\)
‘physical’	\(W - Y^n\)	\(Y^n\)	\(W\)

• eqn=3:

units=	term=’LHS’	term=’RHS’
‘rescaled’	\(\hat{Y} -\left(\frac{m}{M_n}\right)^2\hat{\pi} - \epsilon \left( \frac{M_n}{\Lambda} \right)^{3n-1} \left( \frac{M_{f1}}{M_n} \right)^{n-1} \hat{\nabla}^2 \hat{W}\)	\(\frac{\hat{\rho}}{M_P}\frac{Mn}{M_{f1}}\)
‘physical’	\(Y - m^2 \pi - \frac{\epsilon}{\Lambda^{3n-1}} \nabla^2 W\)	\(\frac{\rho}{M_P}\)

units=	term=1	term=2	term=3	term=4
‘rescaled’	\(\hat{Y}\)	\(-\left(\frac{m}{M_n}\right)^2\hat{\pi}\)	\(- \epsilon \left( \frac{M_n}{\Lambda} \right)^{3n-1} \left( \frac{M_{f1}}{M_n} \right)^{n-1} \hat{\nabla}^2 \hat{W}\)	same as ‘RHS’
‘physical’	\(Y\)	\(-m^2 \pi\)	\(- \frac{\epsilon}{\Lambda^{3n-1}} \nabla^2 W\)	same as ‘RHS’

Equation (3) allows to look at the terms in the equation of motion, Eq. (1) and (2) enforce consistency relations \(Y = '\nabla^2\pi\) and \(W=Y^n\).

Note

In this method, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Parameters

eqn

(integer) choice of equation

term

choice of term: for eqn=1,2, valid terms range are ‘LHS’, 1, 2; for eqn=3, valid terms are ‘LHS’, ‘RHS’ and integers from 1 to 4

norm

‘L2’, ‘linf’ or ‘none’. If ‘L2’ or ‘linf’: compute the \(L_2\) or \(\ell_{\infty}\) norm of the residual; if ‘none’, return the full term over the box - as opposed to its norm.

units

rescaled (default) or physical; choice of units for the output

output_label

if True, output a string with a label for the term (which can be used, e.g. in plot legends)

strong_residual_form(sol, units)

Computes the residual with respect to the strong form of the equations.

The total residual is obtained by summing the residuals of all equations:

\[F = F_1 + F_2 + F_3\]

where, in dimensionless in-code units (units=’rescaled’):

\[ \begin{align}\begin{aligned}& F_1(\hat{\pi},\hat{W},\hat{Y}) = \hat{\nabla}^2 \hat{\pi} - \hat{Y}\\& F_2(\hat{\pi},\hat{W},\hat{Y}) = \hat{W} - \hat{Y}^n\\& F_3(\hat{\pi},\hat{W},\hat{Y}) = \hat{Y} - \left( \frac{m}{M_n} \right)^2 \hat{\pi} - \epsilon \left( \frac{M_n}{\Lambda} \right)^{3n-1} \left(\frac{M_{f1}}{M_n}\right)^{n-1} \hat{\nabla}^2 \hat{W} - \frac{\hat{\rho}}{M_P} \frac{M_n}{M_{f1}}\end{aligned}\end{align} \]

and in physical units (units=’physical’):

\[ \begin{align}\begin{aligned}& F_1(\pi,W,Y) = \nabla^2\pi - Y\\& F_2(\pi,W,Y) = W - Y^n\\& F_3(\pi,W,Y) = Y - m^2 \pi - \frac{\epsilon}{\Lambda^{3n-1}} \nabla^2 W - \frac{\rho}{M_P}\end{aligned}\end{align} \]

Note

In this function, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Note that the weak residual in weak_residual_form() is just the scalar product of the strong residuals by test functions.

Parameters

sol

the solution with respect to which the weak residual is computed.

units

‘rescaled’ (for the rescaled units used inside the code) or ‘physical’, for physical units

weak_residual_form(sol)

Computes the residual with respect to the weak form of the equations:

\[F = F_1 + F_2 + F_3\]

with

\[ \begin{align}\begin{aligned}\begin{split}F_1(\hat{\pi},\hat{W},\hat{Y}) & = - \int\hat{\nabla}\hat{\pi}\hat{\nabla}v_1 \hat{r}^2 d\hat{r} - \int \hat{Y} v_1 \hat{r}^2 d\hat{r} \\\end{split}\\\begin{split}F_2(\hat{\pi},\hat{W},\hat{Y}) & = \int \hat{W} v_2 \hat{r}^2 d\hat{r} - \int \hat{Y}^n v_2 \hat{r}^2 d\hat{r} \\\end{split}\\F_3(\hat{\pi},\hat{W},\hat{Y}) & = \int \hat{Y} v_3 \hat{r}^2 d\hat{r} - \int \left( \frac{m}{M_n} \right)^2 \hat{\pi} v_3 \hat{r}^2 d\hat{r} + \\& + \epsilon \left( \frac{M_n}{\Lambda} \right)^{3n-1} \left(\frac{M_{f1}}{M_n}\right)^{n-1} \int\hat{\nabla} \hat{W} \hat{\nabla} v_3 \hat{r}^2 d\hat{r} - \int \frac{\hat{\rho}}{M_p}\frac{M_n}{M_{f1}} v_3 \hat{r}^2 d\hat{r}\end{aligned}\end{align} \]

The weak residual is employed within solver.Solver.solve() to check convergence - also see solver.Solver.compute_errors().

Parameters

sol

the solution with respect to which the weak residual is computed.

Gravity

class phienics.gravity.PoissonSolver(fem, source, fields=None, Mp=1.0, Mn=None, abs_du_tol=None, rel_du_tol=None, abs_res_tol=None, rel_res_tol=None, max_iter=None, criterion=None, norm_change=None, norm_res=None)

Computes the Newtonian potential \(\Phi_N\) by solving the Poisson equation:

(5)\[\nabla^2 \Phi_N = \frac{\rho}{2 {M_P}^2}\]

where \(\rho\) the source and \(M_P\) the Planck mass.

get_Dirichlet_bc()

Returns the Dirichlet boundary conditions for gravity (Eq. (5)):

\[\Phi_N(\hat{r}_{\rm max}) = 0\]

output_term(term='LHS', norm='none', units='rescaled', output_label=False)

Outputs the left- and right-hand side of the Poisson equation Eq. (5).

The terms can be output in physical or dimensionless in-code units, by choosing either units=’physical’ or units=’rescaled’:

units=	term=’LHS’	term=’RHS’
‘rescaled’	\(\hat{\nabla}^2\Phi_N\)	\(\frac{1}{2} \frac{\hat{\rho}}{M_P} M_n\)
‘physical’	\(\nabla^2 \Phi_N\)	\(\frac{\rho}{2 {M_P}^2}\)

where \(\Phi_N\) is the Newtonian potential, \(\rho\) the source and \(M_P\) the Planck mass.

Note

In this method, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Parameters

term

choice of ‘LHS’ and ‘RHS’ for the left- and right-hand side of the Poisson equation

norm

‘L2’, ‘linf’ or ‘none’. If ‘L2’ or ‘linf’: compute the \(L_2\) or \(\ell_{\infty}\) norm of the residual; if ‘none’, return the full term over the box - as opposed to its norm.

units

rescaled (default) or physical; choice of units for the output

output_label

if True, output a string with a label for the term (which can be used, e.g. in plot legends)

solve()

Overrides the base-class nonlinear solver with a linear solver for the Poisson equation Eq. (5).

strong_residual_form(sol, units)

Computes the residual \(F\) with respect to the strong form of the Poisson equation Eq. (5).

In dimensionless in-code units (units=’rescaled’) it is:

\[F = \hat{\nabla}^2 \Phi_N - \frac{\hat{\rho} M_n}{2 {M_P}^2}\]

and in physical units (units=’physical’):

\[F = \nabla^2 \Phi_N - \frac{\rho}{2 {M_P}^2}\]

Note

In this function, the Laplacian \(\hat{\nabla}^2\) is obtained by projecting \(\frac{\partial^2}{\partial\hat{r}^2} + 2\frac{\partial}{\partial\hat{r}}\). As such, it should not be used with interpolating polynomials of degree less than 2.

Parameters

sol

the solution with respect to which the weak residual is computed.

units

‘rescaled’ (for the rescaled units used inside the code) or ‘physical’, for physical units