Using machine learning to optimise chameleon fifth force experiments

The chameleon is a theorised scalar field that couples to matter and possess a screening mechanism, which weakens observational constraints from experiments performed in regions of higher matter density. One consequence of this screening mechanism is that the force induced by the field is dependent on the shape of the source mass (a property that distinguishes it from gravity). Therefore an optimal shape must exist for which the chameleon force is maximised. Such a shape would allow experiments to improve their sensitivity by simply changing the shape of the source mass. In this work we use a combination of genetic algorithms and the chameleon solving software SELCIE to find shapes that optimise the force at a single point in an idealised experimental environment. We note that the method we used is easily customised, and so could be used to optimise a more realistic experiment involving particle trajectories or the force acting on an extended body. We find the shapes outputted by the genetic algorithm possess common characteristics, such as a preference for smaller source masses, and that the largest fifth forces are produced by small `umbrella'-like shapes with a thickness such that the source is unscreened but the field reaches its minimum inside the source. This remains the optimal shape even as we change the chameleon potential, and the distance from the source, and across a wide range of chameleon parameters. We find that by optimising the shape in this way the fifth force can be increased by 2.45 times when compared to a sphere, centred at the origin, of the same volume and mass.


Introduction
Through the use of a conformal transformation of the metric tensor, some modifications to the Einstein-Hilbert action can be recast as Einstein's general relativity (GR) but with a new scalar degree of freedom that couples to matter [1].Such a scalar field could be the source of the 'dark energy' we observe through the current expansion rate of the universe [2].This field would also mediate a new fundamental force due to its coupling to matter.This 'fifth force' would be observed as deviations in gravitational phenomena from those predicted by GR.However, to date, all experimental results have been in agreement with the predictions of GR, placing strong constraints on the couplings between such a scalar field and matter, in models with a linear equation of motion [3][4][5].Consequently these fields are no longer viable dark energy candidates.This has motivated research into nonlinear scalar field models that possess screening mechanisms.These models have fifth forces that are suppressed in and around regions of relatively high density, resulting in weaker experimental constraints.Some examples include the symmetron, chameleon, and Galileon models [6][7][8].
In this work we will focus on the chameleon, which is a class of models where the field has an effective mass that scales with the local matter-energy density [6].Consequently, in regions of high density, such as close to the Earth, the field's mass is larger and interactions are suppressed.Despite this, numerous experiments, including torsion balance [9], Casimir [10][11][12][13], levitated force sensors [14], atom interferometry [15][16][17][18][19], atomic spectroscopy [20,21], neutron bouncing [22][23][24][25][26], and neutron interferometry [27,28], have been able to constrain the chameleon parameter space as shown in Figure 1.For a comprehensive review of the constraints on the chameleon model see Refs.[29][30][31].  .The current constraints (shaded regions) on the parameter space of the chameleon model with a V (ϕ) = Λ 5 /ϕ potential.M controls the strength of the coupling to matter and Λ the strength of the self-interactions of the scalar field.The black line indicates where Λ is equal to the dark energy value, Λ DE = 2.4 meV.The red line indicates where the chameleon can explain the difference in the measured and theoretical values of the muon g − 2 value [32].LFS stands for Levitated Force Sensor.Figure reproduced with permission from Ref. [13].
The nonlinearity of the chameleon's equations of motion makes calculating analytic solutions very difficult.Examples do exist in the literature but these are only for highly symmetrical source shapes and often involve approximations that limit the solution to specific regions of parameter space [6,33,34].Constraints from experiments and observations will either come from one of these analytic approximations or from numerically solving the equations of motion for the system being studied.The question may then arise: is the experimental system optimal for detecting chameleons?
We know that the approximate analytic field profiles around spherical and cylindrical sources, as expressed in [6] and [34] respectively, have different forms.Consequently, for sources of equal density and radial size, the resulting fifth force induced by the scalar field will differ even at the same radial distance from the source 1 .This means that in an experiment measuring fifth forces, choosing to use a spherical source over a cylinder will improve the experimental sensitivity.This creates an interesting contrast to Newtonian gravity where Gauss' law requires that the force experienced by a test particle is independent of the matter distribution contained in some exterior surface.Since a dependence on the source shape exists for the chameleon [34,[36][37][38], there must exist a shape that maximises the force when taking all other parameters to be fixed.Using such an optimised source shape in current and near future experiments would yield improvements to experimental sensitivity and thereby strengthen the resulting model constraints.The aim of this work is to demonstrate a method by which this optimised source shape can be found.The exact shape will depend on what observable is being optimised and so will depend on the specific experiment.
In this work we will determine the rotationally symmetric shape which maximises the force at a single point a fixed distance from the source, inside a spherical vacuum chamber.We employ several parameterisations to mathematically describe our source shapes.This reduces the problem to finding the set of parameters that maximise the fifth force, although the dimension of this parameter space will inevitably be large if we wish it to encompass complex shapes.Using the software package SELCIE (https://github.com/C-Briddon/SELCIE.git)[35], meshes were generated of the source shape inside a spherical vacuum chamber, with an additional surface for measurements at a fixed distance from the source.2SELCIE was also used to numerically solve for the chameleon field profile of the system and solve for the field gradient magnitude which is proportional to the force.We then iterate over the mesh vertices along the 'measuring surface' and output the maximum value of the field gradient magnitude.The goal then is to find the set of parameters defining the source shape that maximise the field gradient magnitude outputted by this algorithm.To achieve this optimisation we used the DEAP software package (https://github.com/deap/deap.git)to automate the design and running of a 'genetic algorithm' [41].These genetic algorithms are a class of machine learning algorithms that use concepts of Darwinian evolution to improve a population of solutions so that it converges to an extremum in the solution space [42,43].They also do not require the problem being optimised to have an analytic form, and are efficient even for problems with a large solution space dimension.
We will begin this work by reviewing the chameleon model in section 2. The details of the genetic algorithm will be explained in section 3.In section 4 we will define the shape parameterisations used and show the results of our investigation using the genetic algorithm when the measuring distance from the source and chameleon potential are fixed.Sections 5 and 6 will then explore the effects of changing the measuring distance and potential, respectively.We will then end with our conclusions from this work in section 7.
Conventions: In this work calculations are performed in natural units (c = 1 and ℏ = 1), and we use the mostly plus sign convention for the metric, so that the Minkowski metric is η µν = diag(−1, 1, 1, 1).We also use the reduced Planck mass M 2 pl = (8πG N ) −1 , where G N is Newton's gravitational constant.Partial derivatives are expressed as f ,x = ∂f ∂x and f ,xy = ∂ 2 f ∂x∂y .

Chameleon model
The chameleon model [6], in the Einstein frame (with metric g µν and metric determinant g), has the action where R is the Ricci scalar derived in the Einstein frame, and V (ϕ) is the potential of the chameleon scalar field ϕ.The dynamics of other matter species, φ m , are described by the Lagrangian L m (φ m , gµν ), where gµν is the Jordan frame metric, which is related to the Einstein frame metric by the conformal transformation In the action shown in equation (2.1), we have assumed that the scalar field has a universal coupling to matter.Taking this coupling to be constant in ϕ, the conformal factor is A(ϕ) = exp (βϕ/M pl ), where β is the dimensionless universal coupling constant [30].
In this work we assume a field potential of the form where Λ is a mass scale and n is an integer, following Ref.[6].Application of the least action principle to equation (2.1), and taking the matter component to be non-relativistic, the equation of motion of the chameleon field is where □ is the d'Alembert operator, and ρ is the matter-energy density.This equation resembles a Klein-Gordan equation, with the right-hand side being the derivative of the field's effective potential.Taking the left-hand side of equation (2.4) to be zero, we get that the field value that minimises the effective potential is . (2.5) The mass of the field at this minimum is then . (2.6) We see from equation (2.6) that the effective mass of the field (and by extension its Compton wavelength through the relation λ = 1/m ϕ ) depends on ρ, and in such a way that as ρ increases so does the mass.This results in the field becoming suppressed in regions of higher density.This is the chameleon's screening mechanism which allows it to avoid constraints from Earth bound and solar system experiments while having a non-negligible contribution on cosmological scales.
In this investigation we solved for the chameleon field profiles around source masses of arbitrary shapes.To accomplish this, we used the software package SELCIE [35] to numerically solve a static dimensionless version of equation (2.4), for the various shapes, using a finite element approach.To derive this dimensionless form we rescale the density and field values as ρ = ρ/ρ 0 and φ = ϕ/ϕ 0 , where ρ 0 is the vacuum density and ϕ 0 is simply equation (2.5) evaluated at ρ = ρ 0 .We will also work using lengths that are rescaled by the size of our domain, L, and can therefore rewrite the Laplacian as L∇ = ∇.The dimensionless chameleon equation then has the form where the dimensionless constant α is defined as . (2.8) The force experienced by an unscreened test particle, with mass m, due to the chameleon field is Using the rescaling convention above, this formula can be written as where we have also rescaled the force by the magnitude of the gravitational force ⃗ F g which is related to the gravitational acceleration by ⃗ F g = m⃗ a g .Therefore, for fixed model parameters, by maximising the rescaled field gradient magnitude, | ∇ φ|, we will also be maximising the corresponding fifth force.
As shown in Ref. [35], the rescaled Compton wavelength of the field at the minimum of the effective potential and when ρ = ρ 0 is . (2.11) Therefore when α ≫ (n + 1) the Compton wavelength will be much greater than unity.Physically this means the field does not have enough space in the unit sized chamber to reach its vacuum value.Consequently, the rescaled effective potential of the field (see righthand side of equation (2.7)) can be approximated by Veff ≈ − φ−(n+1) inside the vacuum region, since φ ≪ 1 and ρ = 1.If we then apply the substitution φ = α −1/(n+2) φ, the resulting equation of motion for φ will be independent of α inside the vacuum region to leading order.The significance of this result is that by solving for | ∇ φ| for one value of α, we can then determine the fifth force strength for all values where α ≫ (n + 1), by applying an overall rescaling to the solution.This statement remains true while the dense regions of the system are much more dense than the vacuum, because the value of φ inside the dense regions will differ for different values of α, but for sufficiently large density the differences become negligible as φ ≈ 0. Therefore, assuming the above conditions to be satisfied we can express equation (2.10) as where | ∇ φ| (which in the rest of this work will be referred to as the force) depends on n and the profile of ρ but not on Λ, β or L. This equation is also entirely independent of both the source and vacuum density values, so long as there exists a sufficiently large hierarchy between them.We also see from this equation the importance of the domain size on the measured fifth force.Since β and Λ are constants, and φ is independent of L (as long as α ≫ (n + 1)), the total dependence of the domain size, L, on the above expression is L −n/(n+2) .This means smaller systems will yield higher sensitivity to fifth forces.For the case where n = 1 equation (2.12) can be used to express the acceleration acting on the test mass due to the scalar field in physical units as (2.13)

Genetic algorithm
In this work we use the genetic algorithm (GA) software DEAP [41] to optimise the shape of a source mass such that the fifth force it gives rise to at a fixed distance is maximised.GAs are designed to take an initially random population of solutions to a problem and evolve them in order to find the optimum solution [42,44].This is done in three stages; selection, crossover, and mutation.In the first stage, selection, it is determined which individuals will survive to parent the next generation.To make this determination each solution in the population is evaluated using the function we wish to maximise (or minimise), which we will refer to as the fitness function.The fitness of an individual in our population is then the corresponding output of this fitness function, and will determine the likelihood of an individual surviving.As is the case with the later stages (crossover and mutation) there are many methods to determine survivors in our population.For a comparison of the different possible configurations of the GA see appendix A. The method we ultimately chose was a tournament method, where we randomly select N tournsize individuals and return the one with the best fitness (largest value if searching for maximum and smallest value for minimum) [45].This process is repeated N times, where N is the size of our original population.Note that it is possible for individuals to be selected multiple times.This process imitates natural selection and results in a new population with an increased average fitness.It is these survivors that will be the parents of the next generation.The fact that an individual can be selected multiple times is analogous to members of a species that have more offspring, thereby making its characteristics more abundant in the next generation.Reproduction is replicated in the crossover stage, whereby individuals in the population are paired and each pair has some probability, P x , of reproducing.If this occurs then each attribute of the solution pair (e.g.their coefficients) have a p x chance of being swapped.The pair has now become a pair of children that replace the parents in the population.This process combined with the selection process allows individuals with better fitness values to share features with a larger number of other individuals at a faster rate than those with worse fitness.Overall, this will lead the population as a whole to flow towards optimum fitness.
A general problem with optimisation algorithms is that of local maxima.If an individual finds such a solution, it can lead the rest of the population to it without finding the global maximum (e.g., see Ref. [46]).The strategy used by a GA to address this is contained in the mutation stage.This stage is inspired by the mutations in natural organisms required for evolution, and modifies the information present in the current population.In this stage, each individual in the population has some probability, P m , of receiving a mutation.This involves changing some of the attributes of the individual through some process involving randomness.In this work, since our parameters are continuous, each coefficient of the mutated individual has probability p m of receiving a Gaussian perturbation with standard deviation σ.This allows individuals to escape from local maxima while introducing more variation into the population, allowing more of the parameter domain to be explored.
In total, this algorithm has 7 parameters that must be set by the user, {N , N tournsize , P x , P m , p x , p m , σ}.In this work {N , N tournsize , P x , P m , p x } = {105, 5, 0.7, 0.1, 0.5}.The values of {p m , σ} will depend on the number and domain of the coefficients.The process by which we decided on the methods described above, and the values of the parameters, is discussed in Appendix A. cell size was CellSizeMax = 0.01.The exception to this was the boundary between the vacuum and the chamber wall which had CellSizeMin = 5 × 10 −4 and DistMax = 0.5 instead, since this improved the runtime while contributing negligibly to the solutions' accuracy.The thickness of the chamber wall was set to 0.05.Using this mesh, SELCIE is then used to calculate the scalar field profile and the force induced by the field.Finally, the magnitude of the force is evaluated at the mesh vertices that form the measuring surface.The largest of these values is then the output of the fitness function.The optimal shape will correspond to the global maximum of this function on the shape parameter space.For shape spaces with a dimensionality of one or two, it is reasonable to perform an iterative search of the parameters.However, as the dimension of our parameterised shape space increases this method becomes impractical.We will therefore use the GA, as described in section 3, to find the shape parameters that maximise the outputted force value.
In the following sections we will define our shape parameterisations and will probe the corresponding space of shapes to find the shapes that maximise the chameleon force.We will then compare these optimal shapes to determine if there are any consistent features among them and whether these are indicative of a superior class of shapes.

Ellipsoids
We will begin our investigation with ellipsoidal shapes (for earlier work see Ref. [34]).Being the simplest deviation from spherically symmetric sources, these shapes are parameterised by the horizontal and vertical axis radii, r x and r y respectively.Since the shape space is only 2-dimensional, it can be thoroughly investigated without the need of the GA.The boundary of the ellipse in the (x, y) plane is (x, y) = (r x sin(θ), r y cos(θ)), for θ = [0, 2π).For our purposes, however, it will be more useful to parameterise the ellipsoids by their axis-ratio, ϵ = r y /r x , and the volume of the ellipsoid where we have used the fact that our source has rotational symmetry around the y-axis.
For such ellipsoids with a volume V and axis ratio ϵ, their axis radii can be expressed as (r x , r y ) = (ζ, ζϵ), where the scale factor ζ is Some examples of ellipsoids parameterised by (V, ϵ) are shown in Figure 2.
Figure 3 shows the magnitude of the resulting force, at a distance d = 0.05 from the source, for the values ϵ = [0.01,2.0] and V = [10 −5 , 0.1].A note-worthy feature of this plot is that slight deviations from spherical sources appear to always result in an increase in the measured force, the location of which is around the pointed region of the ellipsoid as predicted in Ref. [37].We also note that Figures 3a and 3b show that for fixed ϵ as V decreases the force increases.However, Figures 3c and 3d show that this trend reverses after some critical volume V c (ϵ).This can be seen clearly in Figure 4 where we have plotted the force against volume for a range of ϵ-values.
Figure 5 shows this critical volume V c (ϵ), at which the force is maximised, for a range of ϵ.In this plot we also show the maximal force obtained for the corresponding ellipsoids.We acknowledge that the plot of V c (ϵ) is not entirely smooth.This is because when performing our simulations we use a mesh of finite precision, resulting in a small numerical uncertainty.Although the relative uncertainty on the force is of the order 10 −4 , the finite mesh precision leads to a larger uncertainty on our measurement of V c .Finally, we acknowledge that our plots would indicate that the ideal ellipsoid lies outside our parameter range.This is because of computational limitations of solving the field equation using meshes capable of resolving such small and elongated shapes.In practise, however, a very small ellipsoid will be impractical in realistic experiments.

Legendre polynomial shapes
Legendre polynomials, P n , are a special set of polynomials that are solutions to the differential equation This set of functions has the property showing that it forms a complete basis in the interval [ −1, 1].This means if we define a set of curves in polar coordinates (R(θ), θ) where and a n are the series coefficients, then in the limit N − → ∞ the set will grow to contain every closed curve that is also symmetric around the y-axis.However, we must still consider how we relate these curves to shapes.For example, for some sets of a n , it is possible for R(θ) to take negative values, which can result in the curve intersecting itself.In this work we take positive r < R(θ) to be inside the shape, and to avoid complications will set R(θ) = R min whenever the value of R(θ) given by equation (4.5) is less than R min .To investigate the shape dependence of the force, we keep the volume of the source fixed.We therefore define a mapping (see Appendix B) from one set of coefficients, a n , to a new set where the corresponding Legendre polynomial shape is similar to the original but with a volume V .Some examples of the kinds of shapes contained inside this class are shown in Figure 6.
It can be shown that the Legendre polynomials, in the domain [−1, +1] have a maximum at P (θ = 0) = 1.Constraining the coefficients a n to be positive, this means the maximum   For practical reasons, we will only consider Legendre polynomial shapes that are contained entirely inside the vacuum chamber, e.g.R max < 1 (after constraining the volume).
Using the method described in section 4 we can obtain a force value from the set of coefficients a n that define a Legendre polynomial shape.The number of coefficients used was N = 20; however, after applying the volume constraint the number of dimensions of the shape parameter space is reduced to 19.Because of the large number of dimensions of the shape space we will use the GA discussed in section 3 to find the shape, parameterised by Legendre polynomials, that maximises the force at a distance d = 0.05 from the source evaluated at a single point.For a volume of V = 0.01 the shape outputted by the GA generated a force of | ∇ φ| d = 5.07 (located along the y-axis), which is a 48% increase when compared to a spherical source of the same volume centred at the origin.The profile of the scalar force induced by the source is shown in Figure 7 along with the position where the force is maximised.We repeated our GA simulation for a range of other volumes, and found that as the volume is decreased the force becomes larger, as was the case with the ellipsoids.We also observed that the 'umbrella'-like feature near θ = 3π/4 in Figure 7 appears to be a common feature in each case, with the extra volume being diverted to the upper region, as seen in Figure 8.Since our populations have independently evolved this feature, it suggests this may be the most important feature of the final shape and that it is independent of the total volume.This aspect will be further investigated in section 4.4.

Radially-parameterised shapes
In this class the N points defining the source boundary are given in polar coordinates as (R n , θ n ), where θ n = πn/N and R n is a product of the coefficients b k such that For the case b k = 1 for all k-values, the shape will approximately be a sphere (once the rotational symmetry is taken into account).To ensure the class contains well-behaved shapes we enforce that b k > 0, so radial values are always positive.Also, since each radial value can be seen as the previous value multiplied by a factor b k , by changing the range of values b k can take, we can control the smoothness of the shapes.This second point is important since we are refining the boundary of the shape, and so a large perimeter will require more cells, increasing the overall simulation time and memory required.We again will constrain the volume of the shapes to be a predetermined value V , therefore reducing the dimension of the shape space to N − 1.This was done using the same method as the Legendre polynomial shapes, and is discussed in Appendix B. Some examples of shapes contained in this class are shown in Figure 9.We used the GA to optimise the force, a distance d = 0.05 from the source, of a radiallyparameterised shape defined by N = 50 coefficients ranging between 0.5 and 1.5, and with a target volume of V = 0.01.The shape outputted by the GA had a force of | ∇ φ| d = 5.16 (which is slightly larger than the value obtained for the optimal Legendre polynomial shape), and was again located along the y-axis.The profile of the force generated by this optimised source shape is shown in Figure 10, along with the position of the recorded value.Comparing this shape to the optimal Legendre polynomial shape shown in Figure 7, we see that the GA has outputted a similar shape, with a large lobe around θ = 0 and a protruding spike (umbrella-like shape) in the lower region near θ = 3π/4, despite being generated from two different shape parameterisations.This similarity is not only present for the V = 0.01 case.Comparing the GA results for radially-parameterised shapes of different volumes, we noticed similar shapes to those found when the GA was used with Legendre polynomial shapes for those same volumes, which can be seen by comparing Figures 8 and 11.In section 4.2 we hypothesised that the umbrella-like feature might be the most important to maximising the force, and that the rest of the shape is just there to satisfy the volume constraint.These new results support this hypothesis since we see the same feature independent of the volume and for V = 10 −4 we obtained only this feature without the upper lobe (note that it is not possible to obtain this shape with our Legendre polynomial parameterisation due to the property shown in equation (4.6)).

Is the umbrella the most impactful feature?
In sections 4.2 and 4.3 we found that for both the Legendre polynomial and radiallyparameterised shapes, the GA appeared to converge to similar shapes.Furthermore, shapes obtained from different volume constraints consistently produced a similar feature, as seen in Figures 8 and 11.As we impose rotational symmetry around the y-axis, the feature of interest is an umbrella shape.To test the importance of this umbrella feature we investigated how the measured force responds when other sections of the shapes, obtained in the previous two subsections, are removed.To do this we removed any part of the shape that lies above a cut-off height y c = δy max , where y max is the maximum y-value of the shape and δ ∈ [0, 1].The resulting forces for the Legendre polynomial and radially-parameterised shapes are plotted in Figures 12a and 12b respectively.In each case we see the position where the force is maximised is unchanged while its value increases as more of each shape is removed.This is in agreement with what was observed in section 4.1 (in Figure 4), in that above some critical volume the force increases with decreasing source volume.This further suggests that the umbrella shape, of this particular size, is the optimal shape we have been looking for.We will now test this hypothesis using a new shape class inspired by the umbrella shape in hopes to further improve our force.
This new class consists of three polynomials (one of order N and two of order N + 1) defined over the range x ∈ [0, L u ], where L u is a parameter of the shape.The first polynomial (P 0 ) acts as the baseline of the shape and is defined to have N roots that are contained in the shape coefficients.The magnitudes of the second (P + ) and third (P − ) act as a displacement between P 0 and the boundary of the shape in the positive and negative y-directions respectively.These two polynomials are also defined by N roots contained in the shape coefficients but also have an additional root at x = L u so that the combined curves form a closed loop.For simple polynomials with not too many curves this will lead to umbrellalike shapes.However, this class still has the possibility of producing much more complicated shapes.For further generality the shape is then rotated around the origin anti-clockwise in the (x, y) plane by an angle θ, and then translated by the vector (x 0 , y 0 ).Altogether, the number of parameters defining this class is 3N + 4. In this and future sections we used polynomials of order 5 and 6, so the number of coefficients is 19, each ranging between 0 and 0.3.Unlike the previous shape classes, we place no constraints on the volume of the shape.We do, however, set a minimum positive value on the polynomials P + and P − , to ensure that the curve does not intersect itself.A depiction of how shapes in this class are constructed, along with some examples, is shown in Figure 13.Among these examples one shows a case where the source has been separated into two disconnected pieces (shape filled in black).This occurred because to impose axis-symmetry we only consider parts of the shape for which x > 0, therefore a sufficiently curvy shape can become two disconnected sources.
Using the GA to find the shape in this class that optimises the force we obtained the shape shown in Figure 14.We see the shape outputted by the GA is in agreement with our hypothesis that the umbrella was the optimal shape.Furthermore, thanks to the extra freedom provided by our parameterisation the outputted shape produces the best force value so far at | ∇ φ| d = 5.77, and with a volume of V = 1.26 × 10 −5 .For scale this force is 2.45 times larger than that generated by a sphere, located at the origin, of equal volume and using the same measuring distance.Furthermore, since the volume (and by extension mass) of the source is so small, this will lead to large ratios between the chameleon and gravitational forces.

Effects of varying measuring distance
In this work so far the distance between where we wish to measure, and therefore maximise, the force and the surface of the source has been fixed at d = 0.05.In this section we investigate how varying this distance, d, affects the optimal shape.We will use the polynomial shape class discussed in section 4.4 as it has provided the best results so far.The shapes shown in Figure 15a are the results of the GA using this class, for varying values of d.We see that the umbrella shape remains the optimal shape found in each case, with the point of greatest force (at our fixed distance) being on the y-axis.It can also be seen that the vertical position and length of the shape does depend on the measuring distance, with both increasing with d.
For each shape in Figure 15a we measured the maximum force at a range of distances, the results of which are shown in Figure 15b.In this plot we have also indicated where on each curve the value d is equal to the measuring distance the corresponding shape was optimised for, d opt .As expected each shape has a force that monotonically increases with decreasing distance; however, for d > d opt the force decreases much faster than when d < d opt .Focusing on the points where d = d opt we see that they approximately follow a 1/d relation, and the largest forces for any value of d are from the shapes optimised at that specific measuring distance.This provides further evidence that the optimal shape does depend on the measuring distance in the experiment.This then means that when designing an experiment the minimum measuring distance (due to the experiment's design or from physical effects such as Casimir or Van der Waals forces) must be established first before the optimal source shape can be determined.

Effects of varying n
In the preceding sections we have explored the source shapes that could optimise the chameleon force measured in an experiment, for the potential shown in equation (2.3) with n = 1.However, this is not the only possible choice, and as more of the parameter space for n = 1 is excluded, interest will increase in the n > 1 chameleon models.To investigate how changing the value of n would affect the optimal source shape, we used the GA to find the shape in the polynomial class (see section 4.4) that optimises the force at a separation distance d = 0.05 for n = 2 and n = 3.We note that n can take some negative values that were not investigated in this work, but we expect the findings of this section to hold for these cases as well.The shapes outputted by the GA, and corresponding force profiles, are plotted in Figure 16a along with the n = 1 result discussed in section 4.4, for comparison.We see that the umbrella again appears as the optimal shape in each case, but with an increasing thickness and vertical displacement with increasing n.
We see that the n = 2 shape has a hole in it, and that the n = 3 shape appears to be an umbrella with a 'spike'.Looking at other shapes in the GA population for the n = 2 simulations we find almost identical shapes that generate very similar forces as the optimal but with no hole.We therefore conclude that the existence of the hole contributes negligibly to the force.This is because the size of the hole is smaller than the field's Compton wavelength.Since the Compton wavelength is the inverse of the field's mass, as expressed in equation 2.6, as n increases so will the Compton wavelength and consequently smaller features will not affect the chameleon field profile and not contribute as significantly to the force.We also see this in the n = 3 case where the shape outputted is consistent with the umbrella-like shape but with a spike.Looking at the force magnitude plotted in Figure 16d we see the spike contributes little to the force.The reason the GA outputted this shape can then be considered a consequence of an increasing Compton wavelength leading to a larger degeneracy in the shape to force mapping.Now that we have determined that both the Compton wavelength and thickness of the shapes are increasing with n, we should determine whether the two effects are related by checking how much screening is occurring in each case.In this work we consider the interior of the source to be screened if the field has reached the value that minimises its effective potential in the source, φc (see equation (2.5)), and the field gradient has a small magnitude, | ∇ φ| < 10 −10 .Under this definition we find that each shape shown in Figure 16 is unscreened.However, measuring the field inside the sources, we found that the minimum field value was of the same order of magnitude as φc .This combination of properties make intuitive sense as an unscreened source has all of its mass contributing to the field, and minimising the field's value inside the source allows for larger gradients around it.The reason then for the thickness of the shapes in Figure 16 to increase with n is because as the field's Compton wavelength increases, so does the thickness needed for the field to reach φc .This also helps to explain why the umbrella shape appears to be the optimal shape, since it maximises the amount of unscreened matter around the measuring position while also maintaining the optimal thickness.

Conclusion
In this work we aimed to show that by exploiting the shape dependence of the chameleon field, we could improve experimental sensitivity through the use of an optimised source shape.We also aimed to determine what characteristics were necessary for such a shape.To accomplish this we used the software package SELCIE to numerically solve for the chameleon field profile of an axis-symmetric source mass inside a spherical vacuum chamber, and to measure the resulting force at a fixed distance from the source.The first class of shapes we investigated using this method were ellipsoids centred at the origin.We found that, for all volumes tested, as the ellipsoids deviate from the spherical limit the scalar force increases, with the largest forces being generated by disc-like ellipsoids.Additionally, we found that if the axis-ratio is fixed, then there is a critical volume which maximises the force, resulting in a preference for very small sources.
We investigated more complex shapes, with fixed volumes, using multiple shape parameterisations combined with a genetic algorithm, using the software package DEAP.We found that for each parameterisation used, the genetic algorithm would converge to shapes with a  common feature, even when the volume of the shapes was varied, indicating a solution that is independent of the choice of parameterisation.After further investigation using a new shape parameterisation that did not constrain the volume of the source, we confirmed that the shapes that generate the largest force values are 'umbrella'-like shapes, with the point where the force is maximised being located on the axis of rotation.The force generated by the optimum shape found (for n = 1) was 2.45 times larger than that of a sphere of equal volume located at the origin when the force is measured at the same distance to the source.We also found that the thickness of the umbrella is tuned so that it was both unscreened and had the field reach the value that minimises its effective potential inside the source.This means all the matter sources the field profile whilst also having as large a total field variation as possible, leading to high field gradients.We found that when changing the form of the field's potential, the thickness of the umbrella changed with the field's Compton wavelength, in such a way that this feature was present in each case tested.This characteristic is consistent with the critical volumes found for ellipsoids.When varying the measuring distance d, we found that the umbrella shape remains the optimal choice but its vertical position and length have a dependence on d.Furthermore, we found that the optimal force varies as 1/d.
The umbrella shapes obtained in this work were designed to maximise the fifth force at a single point inside a spherical vacuum chamber.It is not necessarily a given, therefore, that they would be the best shape in an actual experiment, as experiments may use either moving particles, such as in atom interferometry, or an extended body as in torsion balance experiments.However, we believe that the property of the umbrella shapes, whereby the source is unscreened but the field reaches the value that minimises the effective potential inside the source, are generally necessary for optimising a fifth force experiment.Additionally, the method used in this work is general enough and sufficiently customisable that it can be tailored to specific experiments.One would only need to change the value outputted by the fitness function, obtained from the field solution, to whatever observable is used in that particular experiment.This is easily accomplished by extracting the required information from the field profile calculated using SELCIE, for example the total work done on a particle trajectory or scalar forces acting on an extended object.Any experimental constraints such as laser paths, particle trajectories, or different vacuum chamber shapes can also be enforced at the level of mesh generation using shape manipulation tools in SELCIE.These types of changes will likely break rotational symmetry and would therefore require solving the chameleon equations of motion using 3D meshes, which will greatly increase the simulation runtime and will therefore be left to future works.parameters.Before using the GA to optimise our fitness function, as laid out in section 3, we first tested the various possible configurations of the GA to find the most optimal one.This was accomplished using a test fitness function where the fitness value is the non-overlapping cross-sectional area of an input and target shape, both created using the Legendre polynomial parameterisation defined in section 4.2.Therefore, the goal of the GA in these simulations is to minimise this fitness value and return coefficients that are as close to the target ones as possible.
Due to the large number of possible configurations, we decided to treat our choice of selection, crossover and mutation methods as being independent from one another.This greatly reduces the number of configurations that require testing.For each of these configurations the GA was run 10 times.The distribution of the number of generations needed to reach convergence and the corresponding fitness values outputted by the GA are both plotted in Figure 17 for a range of configurations.Additionally, we used this same test to evaluate other in-build algorithms available in DEAP, information on which can be found at https://deap.readthedocs.io/en/master/api/algo.html.Taking both the fitness distribution and number of generations until convergence into account we ultimately found the best choice of those tested was a custom algorithm that uses the methods: cxUniform(indpb=0.5) for crossover (configuration 3 in top plot of Figure 17), mutGaussian(sigma=0.1) for mutation (configuration 1 in middle plot of Figure 17), and selTournament(tournsize=10) for the selection method (configuration 1 in bottom plot of Figure 17).Finally, we tested the effects of changing the input parameters of the chosen methods and found little effect on the distributions of the fitness values or number of generations until convergence.

B Rescaling shapes to constrain volume
In our simulations the continuous space is replaced by a mesh consisting of discrete cells.For a 3D mesh these cells would consist of tetrahedrons with 4 vertices.In this case to measure the volume of our shape we would simply sum the volume of each individual cell that comprises it.For systems with an imposed symmetry, the system can be represented using a 2D mesh consisting of triangular cells made with three vertices.We therefore wish to know how to calculate the volume of the full 3D shape using its 2D intersection and knowledge of the symmetry.
For rotational symmetry we can imagine projecting each of the triangular cells comprising the shape around the axis of symmetry.For cells not touching the axis this will lead to triangular tori centred on the axis, while cells touching the axis will form cone segments.We find that the volume of one of these rotationally extended triangles is where A is the area of the triangle and δ i , for i ∈ {1, 2, 3}, are the distances between the vertices of the triangle and the axis of rotation.To obtain the total volume of the shape, we iterate over every cell comprising the shape and calculate the volume contribution of each cell using equation B.1.In this work, shapes are constructed by joining a series of points together by straight lines to form their boundary.These points are defined on a plane by the radial values R n , for n ∈ [0, N ], and have an angular separation δθ = π/N .Therefore the nth triangle will be defined by the vertices at (R n , nπ/N ), (R n+1 , (n + 1)π/N ), and the origin.Using equation The fitness value is the non-overlapping cross-sectional area of the outputted shape and a target shape with coefficients [0.5, 2, 0.9, 1.3], both of which are created using the Legendre polynomial parameterisation described in section 4.2.In each plot the selection, crossover, and mutation stages not being varied are set as selTournament(tournsize=10), cxUniform(indpb=0.5), and mut-Gaussian(sigma=0.1), respectively.The number of individuals in each population was 100.In the top plot the crossover method is varied with the configuration number, from left to right, corresponding to [cxOnePoint(), cxTwoPoint(), cxUniform(indpb=0.5),cxBlend(alpha=0.1),cxSimulat-edBinary(eta=0.1), cxSimulatedBinary(eta=0.5), cxSimulatedBinaryBounded(eta=0.1),cxSimulat-edBinaryBounded(eta=0.5)].In the middle plot the mutation method is varied and the configurations are [mutGaussian(sigma=0.1), mutGaussian(sigma=0.3), mutGaussian(sigma=0.5), mutGaussian(sigma=0.7), mutGaussian(sigma=0.9), mutShuffleIndexes(), mutPolynomialBounded(eta=0.1),mutPolynomialBounded(eta=0.By performing the rescaling R i → ϵR i , we see from equation (B.2) that the value transforms as V → ϵ 3 V .This means if we first calculate the volume using equation (B.2) and define ϵ = (V / V ) 1/3 , then by rescaling the radial distances by ϵ this will act to fix the volume of the shape to be the target volume V .

Figure 1
Figure 1.The current constraints (shaded regions) on the parameter space of the chameleon model with a V (ϕ) = Λ 5 /ϕ potential.M controls the strength of the coupling to matter and Λ the strength of the self-interactions of the scalar field.The black line indicates where Λ is equal to the dark energy value, Λ DE = 2.4 meV.The red line indicates where the chameleon can explain the difference in the measured and theoretical values of the muon g − 2 value [32].LFS stands for Levitated Force Sensor.Figure reproduced with permission from Ref. [13].

Figure 2 .
Figure 2. Boundary curves of elliptical shapes of varying volume and axis-ratio.The left plot contains ellipses with ϵ = 0.5, while the right contains ellipses with ϵ = 1.5.

Figure 3 .
Figure 3.The measured maximum force (at distance 0.05 from the surface of an ellipsoid) against axis ratio ϵ = r y /r x .The colour map represents the volume of the ellipsoid.Each subplot shows the results for a different order of magnitude in volume, with the largest volumes in the top left and the smallest in the bottom right.

Figure 4 .
Figure 4. Measured maximum force (at distance 0.05 from the surface of the ellipsoid) against source volume, for a range of axis-ratio values.

Figure 5 .
Figure 5. Critical volume (blue) and corresponding force measurement (red) against axis-ratio, ϵ.The vertical gray dashed line at ϵ = 1 indicates the location of the spherical case.

Figure 6 .
Figure 6.Examples of Legendre polynomial shapes defined using N = 20 coefficients and with a volume fixed to V = 0.01.

Figure 7 .
Figure 7.The chameleon force profile of an N = 20 Legendre polynomial source (boundary of which is indicated by the red line) inside a spherical vacuum chamber.The left plot contains the full solution while the right is zoomed in on the central region.In the right plot the red cross on the y-axis below the source, at around y = −0.075,indicates the position where the force (at a distance from the source of d = 0.05) is maximised.The value of this maximal force is | ∇ φ| d = 5.07.

Figure 8 .
Figure 8. Shapes outputted by the GA when finding the optimal Legendre polynomial shape with N = 20 coefficients.Each line indicates a shape with a different volume, as shown in the legend, from V = 10 −4 to V = 10 −2 .The left plot shows the full shapes, while the right is a zoomed in image of the area inside the dashed line.

Figure 9 .
Figure 9. Examples of radially-parameterised shapes defined using N = 50 coefficients and with a volume fixed to V = 0.01.

Figure 10 .
Figure 10.The chameleon force profile of a radially-parameterised source (boundary of which is indicated by the red line) inside a spherical vacuum chamber.The left plot contains the full solution while the right is zoomed in on the central region.The red cross on the y-axis below the source, at around y=-0.075, indicates the position where the force (at a distance from the source of d = 0.05) is maximised.The value of this maximal force is | ∇ φ| d = 5.16.

Figure 11 .
Figure 11.Shapes outputted by the GA when finding the optimal radially-parameterised shape with N = 50 coefficients that range between 0.5 and 1.5.Each line indicates a shape with a different volume, as shown in the legend, from V = 10 −4 to V = 10 −2 .The left plot shows the full shapes, while the right is a zoomed in image of the area inside the dashed line.

Figure 12 .Figure 13 .
Figure 12.The force around Legendre polynomial (a) and radially-parameterised (b) shapes when parts of the upper lobe are removed.Each curve in plot (a) corresponds to one shape shown in Figure 8, while in plot (b) the curves correspond to the shapes in Figure 11.The volume of the unmodified shapes is shown in the legend.

Figure 14 .
Figure 14.The chameleon force profile of a source whose shape belongs to the polynomial shape class where N = 5 (the boundary of the shape is indicated by the red line), inside a spherical vacuum chamber.The left plot contains the full solution while the right is a zoomed in version showing the region immediately around the shape.The red cross on the y-axis below the source, at around y = 0.09, indicates the position of where the force (at a distance from the source of d = 0.05) is maximised.

Figure 15 .
Figure 15.Each line in plot (a) depicts the outline of a cross section of a polynomial shape (as described in section 4.4) outputted by the GA when optimising the force at different measuring distances, d opt .For each of these shapes we then measured the maximum force at a range of distances, d, from the surface of the shape and plotted the resulting curves in plot (b).The black circles indicate when d = d opt , and the black line shows the curve | ∇ φ| = 0.295/d.

Figure 16 .
Figure 16.Polynomial shapes outputted by the GA for fixed measuring distance d = 0.05, but for varying n.Figure (a) shows the shapes compared to one another.Figures (b), (c), and (d) each show the chameleon force profiles sourced by the shapes in (a), for n = 1, n = 2, and n = 3 respectively.The red lines indicate the boundary of the shape.
Figure 16.Polynomial shapes outputted by the GA for fixed measuring distance d = 0.05, but for varying n.Figure (a) shows the shapes compared to one another.Figures (b), (c), and (d) each show the chameleon force profiles sourced by the shapes in (a), for n = 1, n = 2, and n = 3 respectively.The red lines indicate the boundary of the shape.