Capacitance Calculation for a Parallel-Plate Capacitor: Maxwell Matrix vs. Theoretical Value

[Takashima1996]

Hello, and thank you for developing such a useful library for radiation detectors.

I’m currently testing it by developing a simple program to simulate a capacitor composed of two rectangular electrodes (each 2 x 2 x 0.1 cm) placed in a vacuum, with a 0.99 cm gap between their nearest sides (details in the last part). The theoretical capacitance of this setup is approximately 0.36 pF, calculated using this online calculator.

However, when I calculate the capacitance matrix using SolidStateDetector (see the attached script) and then derive the capacitance using Eq. (34) from this paper, I get around 0.83 pF.

function run_calculation(geom_filepath)
    T   = Float32
    sim = Simulation{T}(geom_filepath)

    apply_initial_state!(sim, ElectricPotential)
    calculate_electric_potential!( sim, refinement_limits = [0.1, 0.05])
    calculate_electric_field!(sim)

    for contact in sim.detector.contacts
        calculate_weighting_potential!(sim, contact.id, 
            refinement_limits = [0.2, 0.1, 0.05], verbose = false)
    end
    return sim
end

function cal_capacitance(C_mat)
    C_cap = (C_mat[1,1]*C_mat[2,2] - C_mat[1,2]*C_mat[2,1]) / (C_mat[1,1] + C_mat[2,1] + C_mat[1,2] + C_mat[2,2])
    return C_cap
end

#Main
geom_filename   = "geometry.yaml"
results, runtime = @timed run_calculation(geom_filename)
cap_matrix       = calculate_capacitance_matrix(results)
cal_capacitance(cap_matrix)

Could anyone help me understand why there’s such a discrepancy between these two capacitance values?

Contents of the geometry yaml file are as follows:

name: Public Capacitor
units:
  length: cm
  angle: deg
  potential: V
  temperature: K
grid:
  coordinates: cartesian
  axes:
    x:
      from: -5
      to: 5
      boundaries:
        left: inf
        right: inf
    y:
      from: -5
      to: 5
      boundaries:
        left: inf
        right: inf
    z:
      from: -5
      to: 5
      boudaries:
        left: inf
        right: inf

medium: vacuum
detectors:
  - semiconductor:
      material: vacuum
      temperature: 87
      geometry:
          box:
            widths: [5.0, 5.0, 0.99] 
            origin: [0, 0, 0]

    contacts:
      - name: "anode"
        id: 1
        potential: 0
        material: copper
        geometry:
          box:
            widths: [2.0, 2.0, 0.01]
            origin: [0, 0, 0.5]
      - name: "cathode"
        id: 2
        potential: -100
        material: copper
        geometry:
          box:
            widths: [2.0, 2.0, 0.01]
            origin: [0, 0, -0.5]

[fhagemann]

Hi, first of all thanks for your interest and for reaching out with this question.
This question is indeed a tricky one, but let me try to answer it.

The mutual capacitance is calculated as given in Eq. (8) of the paper you were quoting and as explained in our documentation:

$$c_{ij} = \epsilon_0 \int_{World} \nabla \Phi_i^w(\vec{r}) ϵ_r(\vec{r}) \nabla \Phi_j^w(\vec{r}) d\vec{r}$$

1. calculate the gradients of the weighting potentials
2. calculate the scalar products of these gradients
3. multiply with ϵ0 * ϵr
4. integrate over the whole world volume

Now, the size of the world chosen for the simulation might influence the result.
In your configuration file, the weighting potentials look like this:

plot(sim.weighting_potentials[1], x = 0, contours_equal_potential = true, levels = 0:0.1:1, linecolor = :white)
plot!(sim.detector, st = :slice, x = 0, label = "", xlims = (-0.05,0.05), ylims = (-0.05,0.05))

Contact 1

Contact 2

Here, not only the volume between the plates contribute to the total capacitance but also the non-vanishing weighting potential gradients, especially left and right of the 2x2x0.99 cm^3 volume enclosed by the contact plates.

The result of calculate_capacitance_matrix(sim) for this setup is:

2×2 Matrix{Quantity{Float32, 𝐈^2 𝐓^4 𝐋^-2 𝐌^-1, Unitful.FreeUnits{(pF,), 𝐈^2 𝐓^4 𝐋^-2 𝐌^-1, nothing}}}:
   1.13642 pF  -0.489627 pF
 -0.489626 pF    1.13615 pF

For your simulations, you might wanna restrict your world to this smaller 2x2x1cm^3 volume (including contacts) by adjusting the grid section in your config file

grid:
  coordinates: cartesian
  axes:
    x:
      from: -1
      to: 1
      boundaries:
        left: reflecting
        right: reflecting
    y:
      from: -1
      to: 1
      boundaries:
        left: reflecting
        right: reflecting
    z:
      from: -0.5
      to: 0.5
      boundaries: # beware, there was a typo in your config file
        left: fixed
        right: fixed

In this case, the weighting potential is only calculated in a part of the previous world

plot(sim.weighting_potentials[2], x = 0, contours_equal_potential = true, levels = 0:0.1:1, linecolor = :white)
plot!(sim.detector, st = :slice, x = 0, label = "", xlims = (-0.01,0.01), ylims = (-0.005,0.005))

Contact 1

Contact 2

with a resulting calculate_capacitance_matrix(sim) of

2×2 Matrix{Quantity{Float32, 𝐈^2 𝐓^4 𝐋^-2 𝐌^-1, Unitful.FreeUnits{(pF,), 𝐈^2 𝐓^4 𝐋^-2 𝐌^-1, nothing}}}:
  0.357745 pF  -0.357745 pF
 -0.357745 pF   0.357745 pF

which agrees perfectly with the quoted value. Small drawback: your method cal_capacitance returns NaN32 pF because of division by zero here.. (the absolute values of all elements are equal..)

In general, it should be noted that the surroundings of the radiation detector can significantly influence its capacitance.
However, if the medium between the plates is some semiconductor (e.g. germanium with ϵr = 16), the contribution from inside the detector will be weighted more heavily than the contribution from the outside.

[Takashima1996]

Thank you for the detailed explanation.

It seems that the discrepancy is indeed due to edge effects from the electrodes. When the distance between the electrodes is set to a much smaller value, the calculated capacitance closely matches that of an ideal parallel-plate capacitor with the same electrode geometry.

I appreciate your help once again.